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ABSTRACT 

The  problem  of  accurately  replicating  the  parameters  which 
define  a  given  system  for  the  purposes  of  implementing  modern 
control  strategies  is  important.   Using  an  Autoregressive- 
Moving  Average  (ARMA)  representation  for  the  unknown  system, 
a  model  is  identified  by  processing  input/output  data  to  esti- 
mate the  coefficients  associated  with  the  ARMA  equation.   Iden- 
tification of  unknown  system  parameters  using  Kalman  filtering 
methods  was  accomplished  by  augmenting  the  state  vector.   In 
this  thesis  the  Kalman  filter  is  formulated  so  that  parameters 
can  be  identified  explicitly.   We  call  this  approach  the 
Adaptive  Kalman  Identifier  (AKI). 

It  is  shown  that  the  Adaptive  least  mean  square  (LMS)/ 
Adaptive  Recursive  LMS  and  Adaptive  Lattice  filters  are  special 
suboptimal  cases  of  the  AKI .   The  convergence  and  modeling 
properties  are  compared  with  those  of  the  AKI  by  simulation 
using  various  types  of  data. 

With  minor  modifications,  the  AKI  algorithm  was  used  to 
identify  the  linear  and  non-linear  ARMA  models  of  the  phase 
locked  loop  (PLL) .   A  discrete  PLL  using  a  forward  Euler  inte- 
gration scheme  was  used  as  a  source  of  non-linear  data.   The 
AKI  technique  appears  to  enable  one  to  discern  when  a  potential 
non-linear  system  enters  into  its  non-linear  mode  of  operation. 
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I.   INTRODUCTION 

The  accomplishments  in  the  area  of  microprocessor  tech- 
nology in  the  last  decade  have  made  a  noticable  impact  in  the 
area  of  modern  control.   The  desire  to  implement  modern  con- 
trol theories  taking  advantage  of  these  advancements,  has 
made  it  imperative  that  the  engineer  attempt  to  mathemat- 
ically replicate  the  system  he  ultimately  desires  to  control. 
Efforts  in  this  regard  have,  hence,  generated  a  growing 
interest  in  the  area  of  system  identification  [Ref.  1,2,3,4]. 
By  implication  this  thesis  concerns  itself  with  discrete/ 
digital  signal  processing. 

A.   BACKGROUND 

System  identification  or  modeling  can  be  accomplished  by 
innovative  application  of  existing  techniques  which  were 
generally  considered  filtering  or  state  estimation  methods 
[Ref.  5].   Previous  research  efforts  have,  for  obvious  reasons, 
focused  on  linear  modeling;  however,  there  is  a  rising  interest 
in  non-linear  modeling  methods  [Ref.  6,7,8]  . 

An  attractive  form  for  the  representation  of  an  unknown 
system  is  the  Autoregressive-Moving  Average  (ARMA)  equation, 

00  CO 

y(k)   =    I      a  u(k-i)  -  I      b.y(k-i)  (1.1) 

i=0   ^         i=l   ^ 

which   states    that   the   present  output,    y(k) ,    is    a   linear   com- 
bination of  past  outputs,    y(k-i) ,    and   of  past   and  present 
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inputs,    u(k).      Its    attractiveness    lies    in   its    linear   charac- 
ter  and  easily    imp lemen table    structure   using   microprocessor/ 
computer    algorithms.       Its    relationship    to    the    state    space 
representation   of    a   linear   system  has    already   been  established 
[Ref.    9,10]    making   it   germane    to   consider    that   a    system    is 
identifiable    by    an   equation   of    the    form    (1.1).       The    system 
identification   problem    thus   entails    identifying    the    coeffi- 
cients   a .    and  b . . 
1  1 

There    are    several  methods    for   computing   the    coefficients, 
a.    and  b- ,    however,    it   is    not   the    intent   of    this   brief   intro- 
duction  to   attempt   to   develop   even   the   majority   of    them. 
Nevertheless,    it   is   practical    to  present   a   referenced  history 
of   those  methods   encountered   in   this    thesis. 

Adaptive    algorithms    for   the   purpose    of  estimating    the 
coefficients    of    (1.1)    have    always   been  of   interest.      Widrow, 
using   a   least  means   square   error   criterion   and   implementing 
the  method  of   steepest  descent,    developed   an   adaptive    algorithm 
which   estimated    the    coefficients    of    the   moving    average    process 
associated  with   equation    (1.1)     [Ref.    11] .       That    is,    the   moving 
average   model. 


y(k)       =      a^u(k)    +aj^u(k-l)   +a'u(k-2)    +    ...  (1.2) 

where , 


^0      "      ^0 


^i    =     ^^1  -  ^0^1^ 
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^2      =       ^^2   "    ^0^2^    -   ^l^^l    -    ^0^1^ 


was    identified.       The    LMS    theory   has    since    been  extended    to 
include    block   LMS    filtering   methods     [Ref.    12,13,14]    using 
various    search    techniques    [Ref.    15,4(Chapter    5)].      These 
methods   have   enjoyed  much  popularity    in   the    area  of   Linear 
Prediction   and   digital    speech   processing    [Ref.    16,17]. 

The    shortcoming   of    representing    equation    (1.1)    by    its 
equivalent  moving   average  model    (1.2)    lies    in    the   practical 
aspect  of    its    implementation.      That   is,    the    infinite    series 
represented  by  equation    (1.2)    must  be    tiruncated   at  some   point 
resulting   in   an   approximation  which  may   not   adequately   repre- 
sent equation    (1.1) .      Hence,    efforts    to   adaptively   estimate 
the    coefficients    for   both    the    autoregressive    (b.)    and  moving 
average    (a.)    processes    continued    [Ref.    18,19,20]    with   various 
degrees    of   success.      Feintuch's   Adaptive    Regressive    LMS   pro- 
cedure   [Ref.    18]    is    still   a   controversial    issue    [Ref.    19,20]. 

From  yet   a  different  direction,    Anderson   and  Moore   suggested 
that   the    Kalman    filtering    algorithms    can   be    used   as    a  means 
to   identify    the    a.    and   b.    coefficients      of   equation    (1.1) 
[Ref.    21,    pp.    50-52].       This    thesis    exploits    this    application. 


Anderson   and  Moore    in   fact   formulate    the    technique    to 
compute    the    coefficients    a.,b.    i   =   1,2,     ...    p  where 
p   =  n+m.  "'"      ■"■ 
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B.   INVESTIGATIONS  AND  CONTRIBUTIONS 

This  thesis  formulates  an  adaptive  application  of  the  Kal- 
man  filter  to  identify  the  coefficients  of  the  general  ARMA 
equation  (1.1).   It  is  shown  that  the  LMS  adaptive  filter 
[Ref.  11]  and  the  Adaptive  Recursive  U4S  filter  [Ref.  18] 
are  (1)  special  cases  of  the  Adaptive  ARMA  Kalman  identifier 
and  (2)  sub-optimal  with  respect  to  the  underlying  least 
means  square  (LMS)  error  criterion  upon  which  the  LMS  adap- 
tive filter,  the  adaptive  recursive  filter  and  the  Kalman 
filter  are  based.   It  is  also  demonstrated  that  the  adaptive 
Kalman  identifier  has  excellent  convergence  properties . 

By  comparison,  it  is  noted  that  the  Kalman  algorithm 
accounts  for  measurement  noise  where  the  LMS  algorithms  do 
not,  a  heretofore  unapproached  problem  by  LMS  algorithms. 
The  results  indicate  that  the  suggested  modification  by 
Griffiths  [Ref.  22]  of  the  convergence  factor,  k  ,  for  the 
LMS  adaptive  filter  is  justified. 

The  application  of  the  Kalman  filter  algorithm  is  ex- 
tended to  identification  of  the  coefficients  associated  with 
a  special  case  of  the  generalized  non-linear  ARMA  model  [Ref. 
8] .   The  results  of  the  non-linear  simulations  suggest  a  tech- 
nique for  determining  when  a  potential  non-linear  system 
enters  its  non-linear  operating  regime.   Such  a  technique  can 
prove  valuable  when  on-line  performance  evaluation  of  a  known 
non-linear  system  is  required. 

Lastly,  the  connection  between  the  Kalman  filter  algorithm 
and  the  lattice  filter  algorithm  [Ref.  6]  is  made.   An  example 
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which   demonstrates    and    compares    the    performance    of   both 
algorithms    is   given. 

C.       ORGANIZATION 

Chapter   II   presents    and  exploits    the    Kalman   filter  equa- 
tions  emphasizing    its    connections   with    the   Yule-Walker   and 
the    discrete   Weiner-Hopf   equations.       The    theory    is    further 
developed    to    investigate    the    general   ARMA   case.       Chapter    III 
pursues    the    theoretical   comparisons    between    the    discrete 
Weiner-Hopf   equation  upon  which    the    adaptive    LMS    filter    is 
based   and    the   MA  form  of    the    Kalman    filter.      The    theoretical 
comparison   between   the   Adaptive   Recursive    LMS    filter   and   the 
general   ARMA   form   of   the    Kalman    filter    is   made    in   Chapter    IV. 
The   software   methods   by  which    linear   and  non-linear   synthetic 
data   is    generated   are   discussed   in  Chapter  V.      A  short  user's 
description   of   the   data  processing   programs    and   the   options 
provided   is    given    in   Chapter  VI.      A  non-linear   application 
for   the    identification   of    the    parameters    of    a   non-linear   plant 
is   presented   in   Chapter  VII    followed  by    a   discussion   and   pre- 
sentation of   the    results    in  Chapter  VIII. 
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II.   FORMULATION  OF  THE  PROBLEM 

Identifying  the  coefficients  of  the  Autoregressive- 
Moving  Average  (ARMA)  equation, 

y(k)  +  b,y(k-l)  +  ...  +  b  y(k:-n)   =   a„u(k)  +  a,u(k-l) 
■'1  n  u        1 

+  ...  +  a  u (k-m)    (2.1a) 
m 

m  n 

y(k)   =   I   a.u(k-i)  -  I   b.y(k-i)  (2.1b) 

i=0  ^         i=l  ^ 

can  be  formulated  as  an  adaptive  Kalman  identification  problem. 
That  is,  instead  of  using  the  well  known  Kalman  Filter  equa- 
tions [Ref.  23,24]  to  recursively  estimate  the  states  of  a 
system,  one  can  utilize  the  Kalman  filter  equations  to  adap- 
tively  estimate  the  coefficients  of  either  an  Autoregressive 
(AR) ,  Moving  Average  (MA) ,  or  Autoregressive-Moving  Average 
(ARMA)  process  by  proper  definition  of  the  quantities  involved. 
We  call  this  the  Adaptive  Kalman  Identifier  for  obvious 
reasons . 

A.   THE  DETERMINISTIC  CASE 

Consider  for  the  moment  that  the  a.  and  b.  are  constant. 

1      1 

Then  by  collecting  a  sufficient  number  of  measurements  of  the 
input  u(k)  and  the  output  y(k) ,  one  can  readily  solve  a  set 
of  linear  equations  for  the  a.  and  b^.  coefficients.   The 
"sufficient  number"  that  is  needed  is  n+m+1  which  define  the 
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n+m+l  linear  equations  which  solve  the  n+m+1  unknown  coeffi- 
cients.  The  matrix  equation  to  be  solved  takes  the  form: 


y(0) 
yd) 


y(m) 
y(m+l) 


y(m+n) 


u(0) 
u(l) 


0 
u(0) 


...  0 


u(m)    u(m-l) 
u(rrH-l) 


u(0) 
u(l) 


u(m+n) 


0 
y(0) 


y(m-l)  y(m-2)    y(0)   0 
y(m)    y(m-l)    y(l)   0 


u(m)  I  y(n-l)  y(n-2) 


y(0) 


^0 


m 


(2.2) 


y   =   T 


(2.3) 


It  is  interesting  to  note  that  the  (n+m+1)  x  (n+m+1)  matrix, 
T,  is  block  symmetric.   The  solution  to  (2.3)  is  readily 
apparent,  namely. 


=  t""^^ 


(2.4) 


where   T        is    the    inverse   matrix   of   T.      Since   we    have    assumed 
that   the   elements    of   T   are   perfect  noiseless  measurements, 
and   that  we    somehow  know   a-priori    the    number  of   unknown 
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coefficients,  then  T  is  of  full  rank,  namely  rank  (T)  = 
n+m+1,  and  the  inverse  of  T  exists. 

B.   THE  NONDETERMINISTIC  CASE 

Rarely  can  the  coefficients  be  modeled  so  ideally.   A 
more  prudent  and  realistic  model  admits  that  the  ARMA  equa- 
tion coefficients  are  subject  to  random  perturbations.   Further, 
it  can  be  said  that  in  general  measurement  devices  introduce 
noise  into  the  measurement  data.   Hence,  the  measurement 
model  should  reflect  this  fact.   Developing  the  Adaptive  Kalman 
Identifier  along  the  same  lines  as  [Ref.  21]  we  let. 


a.  (k+1)   =   a.  (k)  +  w.  (k)     i   =   0,1,2,  ...  m    (2.5a) 

b.(k+l)   =   b.(k)  +w.(k)     j   =   1,2,  ...  n      {2.5b) 

where  the  {w. (k),w.(k)}  are  samples  from  a  zero  mean,  white, 
gaussian  random  process.   Additionally,  we  assume  that  the 
noise  sources  are  uncorrelated. 


E[w   (k)w   (k)]   =   0   for   r  7^  s  (2.6a) 

a     a 
r     s 

E[w,   (k)w,   (k)]   =   0   for   r  7^  s  (2.6b) 

b     b 
r      s 

E[w   (k)w,   (k)]   =   0   for  all  r,s  (2.6c) 

a     b 
r     s 

Similarly,  allowing  for  noisy  measurement  devices,  the 
measurement  equation  (2.1b),  is  modelled  as,' 


m  n 

y(k)   =    I      a.u(k-i)  -  I      b.y(k-i)  +  v(k)        (2.7) 
i=0   ^         i=l   ^ 
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where  the  {v(k)}  are  samples  from  a  zero  mean,  white,  gaussian 
random  process.   It  is  also  assumed  that. 


E[w. (k)  v(k) ]   =   0   for  all  i 


E[w.  (k)  v(k)  ]   =   0   for  all  j 


(2.8a) 
(2.8b) 


That  is,  it  is  assumed  that  the  noise  perturbations  associated 
with  one  coefficient  are  independent  of  the  noise  perturba- 
tions associated  with  any  other  coefficient,  and  that  the 
measurement  noise  and  the  coefficient  perturbations  are 
independent. 

At  this  point  it  is  ncessary  to  use  judgment  and  experi- 
ence and  utilize  all  the  information  known  about  the  physical 
system  which  equation  (2.7)  represents  to  assign  variances  for 
the  random  processes  {w.(k)},  {w.(k)}  and  {v(k)}.   Let 


Q   =  E 


w^(k) 


w  .  (K) 


[W^(k)   Wj(k)]   > 


Q   = 


2  2 

where,  Q,  =  diag  (a   ) ,  Q^  =  diag  (a   )  and, 
1  w .     z  w . 

1  3 


R   =   E[v(k)  v(k)] 


(2.9a) 


(2.9b) 


(2.9c) 


and  E[x]  is  the  expected  or  average  value  of  x. 
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Instead  of  using  x's  denoting  the  commonly  used  notation 
for  state  variables,  we  retain  the  flavor  of  the  problem  by 
using  the  a.'s  and  b. 's  as  the  "states"  of  the  Adaptive  Kalman 
Identifier.   It  is  hoped  that  a  more  meaningful  understanding 
of  the  Adaptive  Kalman  Identifier  may  thus  be  gained.   There- 
fore, define  the  state  vector, 


ao(k) 
a^(k) 


m 


b2(k) 


(2.10) 


b^(k) 


Combining  equations  (2.5)  and  (2.10),  we  have  the  discrete, 
first  order  Gauss-Markov  process  [Ref.  3], 


a(k+l) 
b(k+l) 


a(k) 
b(k) 


+  w(k) 


(2.11) 


whe  re , 


w(k)   = 


w . 

— 1 

(k)" 

w  . 

(k) 

(2.12) 
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Note  that  w(k)  is  an  n+m+1  vector  whose  elements  are  a  con- 
catenation of  the  noise  sequences  associated  with  the  a.  and 
b..   To  complete  the  Adaptive  Kalman  Identifier  formulation 
we  define  the  measurement  vector. 


H(k)   =   [u(k)  u(k-l)  ...  u(k-m)   -y(k-l)  ...  -y (k-n) ]   (2.13) 


and  its    associated  measurement  equation, 


y(k)       =      H(k) 


a(k) 
b(k) 


+  v(k) 


(2.14) 


Then  the  solution  to  the  Adaptive  Kalman  Identifier  problem 
[Ref.  21]  is. 


a(k+llk) 


b(k+l!k) 


[I  -K(k)H(k)  ] 


a(k|k-l) 


b(k  ik-l) 


+  K(k)y(k)     (2.15a) 


K(k)   =   P(klk(l)  H'^(k)  [H(k)  P(k!k-1)  H'^(k)+R]"-^ 


(2.15b) 


P(k+llk)   =   P(klk-l)  -  K(k)H(k)P(k|k-l)  +Q. 


(2.15c) 


Equations  (2.15)  are  initialized  by  assuming  an  initial  value 
for  the  coefficients  (2.10)  and  assigning  to  our  assumption 
a  measure  of  our  confidence  in  the  initial  guess.   That  is, 
we  pick  the  values , 


a(0|-l) 
b(Ol-l) 


and  P(0  -1) 


(2.16) 
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where  P(0|-1)  is  defined  as  the  error  covariance  of  the 
coefficients , 


P(klk-l)   =   E 


a(k) 
b(k) 


a(k  ik) 


b(k|k) 


a(k) 


b(k) 


a(k   k) 

T  ' 

b(klk) 

. 

(2.17) 


for  k  =  0.   Equation  (2.15b)  is  generally  referred  to  as  the 
Kalman  gain.   Two  special  cases  are  of  interest:   case  (1); 
a.  =  0  for  all  i,  and  case  (2) ;  b.  =  0  for  all  i. 
1.   Autoregressive  Form  (a.  =  0) 


For  case  1,  equation  (2.7)  takes  the  form, 


y(k) 


n 

y   b.y(k-i)  +  v(k) 
i  =  l   ^ 


(2.18) 


This  is  a  recursive  equation  stating  that  the  present  output 
is  a  linear  combination  of  past  outputs  corrupted  with  addi- 
tive gaussian  noise,  v(k).  Equation  (2.18)  is  more  formally 
recognized  as  an  autoregressive  process  of  order  n  [Ref.  25, 
26] .  Writing  the  Kalman  solution  for  the  coefficients  of  the 
AR  process, 


b(k+l  Ik) 


[I  -K(k)K(k)  ]b(k|k-l)  +K(k)y(k) 


(2.19a) 


K(k)   =   P(klk-l)  H'^(k)  [H(k)P(k  Ik-l)  H'^(k)+R]"^     (2.19b) 
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P(k+l|k)   =   P(klk-l)  -  K(k)H(k)P(k!k-l)  +Q.         (2.19c) 

Alternate  forms  of  the  Kalman  equations  given  by  Maybeck 
[Ref.  27]  are, 

b(k+l|k)   =   P(k|k-1)  P(k|k-1)"-^  b(klk-l) 

+  P(k+l|k)  H'^(k)  R"-^y(k)  (2.20a) 


P(k+llk)   =   (P(k  Ik-l)  )"^  +H'^(k)  r""'-  H  (k)  ""'■ .  (2.20b) 


We  can  model  the  fact  that  we  have  no  a-priori  knowledge 
about  the  initial  values  of  the  coefficients  by  letting, 

[P(0|-1)]~-^   :::   0.  (2.21) 

Further,  if  we  remain  ignorant  and  totally  doubt  our  previous 
estimate,  then  P(k|k-1)  is  modeled  as. 


[P(k|k-1)  ]"-^  z      0.  (2.22) 


Equations  (2.20)  reduce  to. 


b(k+l|k)   =   [H(k)  r"-^  H'^(k)]"^  H(k)  R"-'-Y(k),      (2.23) 


the  weighted  least  squares  estimate,  previously  encountered 
by  Swerling  [Ref.  28,29,30],  of  the  coefficients,  b.   Carry- 
ing the  analysis  further  and  letting  R   =  1 ,  we  have  the 
Penrose  pseudo-inverse  solution  [Ref.  31,32], 

b(k+l|k)   =   [H(k)  H'^(k)  ]"-^  H(k)y(k)  (2.24) 
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Swerling  [Ref.  28]  has  shown  that  if  the  weighting  matrix 
in  (2.23)  is  not  the  inverse  of  the  covariance  matrix  of 
the  measurement  errors ,  then  the  accuracy  of  the  estimated 
coefficients  (2.24)  will  be  degraded. 

The  aforementioned  notwithstanding,  we  press  further 

T 
into  the  analysis  of  equation  (2.24).   The  product  H(k)H  (k) 

can  be  written  as, 

y(k-l) 
y(k-2) 


H(k)H"(k) 


y(k-n) 


[y(k-l)  y(k-2)  ...  y(k-n)] 


(2.25) 


H(k)r(k)  = 


y  (k-1) 


y(k-l)y(k-2)  ...  y(k-l)y(k>n) 


y(k-2)y(k-l)  y"(k-2) 


y(k-n)y(k-l) 


...  y(k-n)y(k-n) 


(2.26) 


As  we  let  k  -*  0°  the  expected  value  of  (2.26)  becomes  the 
autocorrelation  matrix. 


Edim  H(k)  H  (k)  }   =   R   (k) 


(2.27) 


where , 
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R   (k)   = 

yy 


R   (0)  R   (1) 

yy  yy 

R   (-1)  R   (0) 

yy  yy 


R   (n) 

yy 


R   (-n) 

yy 


. . .   R   (0) 

yy 


similarly,  the  product  H(k)y(k)  as  k  -^  =°  becomes 


R   (-1) 

yy 


r   (k: 

yy 


R   (-2) 

yy 


(2.28) 


(2.29) 


R   (-n) 

yy 


The  steady  state  solution  for  the  estimate  of  the  coeffi- 
cients  b(k+l|k)  is, 


b   (k+llk)   =   R   "■^(k)r   (k)  . 

— ss    '      yy     yy 


(2.30) 


Equation  (2.30)  is  one  of  the  starting  points  from  which 
Perry  [Ref.  6]  develops  his  Lattice  modeling  algorithms. 
Equation  (2.30)  can  also  be  recognized  as  the  solution  to 
the  Yule-Walker  or  Normal  equations  [Ref.  26]. 
2.   Moving  Average  Form  (b.  =  0) 


Referring  once  again  to  equation  (2.7)  and  setting 
the  coefficients  b.  =0  for  all  i  we  have  case  (2) , 


m 


y(k)   =    I      u(k-i)  +  v(k) 
i=0 


(2.31) 
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since  v(k:)  is  by  definition  noise  associated  with  the 
measurement,  we  can  combine  y(k)  and  v(k)  such  that. 


m 
z(k)   =   y(k)  -  v(k)   =    I      a.u(k-i).  (2.32) 

i=0   ^ 


Equation  (2.32)  simply  states  that  the  present  measurement 
is  a  linear  combination  of  past  and  present  inputs  or  by 
definition,  a  Moving  Average  (MA)  process  [Ref.  25,33].   The 
Adaptive  Kalman  Identifier  estimate  for  the  coefficients  of 
the  MA  process  is  as  before, 

a(k+llk)   =   [I-K(k)H(k)]  a(k+llk)  +K(k)z(k)         (2.33a) 


K(k)   =   P(klk-l)  H^(k)  [H(k)P(k  Ik-l)  H*^  (k)  +  r]  ""'■      (2.33b) 


P(k+l|k)   =   P(k|k-1)  -K(k)H(k)P(klk-l)  +Q.  (2.33c) 

Using  the  alternate  forms  of  the  above  equations  we  arrive 
at. 


a(k+l|k)   =   [P(k+l!k)  [P(k  |k-l)  ]"-^]a,(klk-l) 

+  [P(k+l!k)  H^(k)  R"-^]z(k)  (2.34a; 


P(k+llk)   =   [  [P(k|k-1)  J"-"-  +H^(k)  r"-*-  H(k)  ]"^.         (2.34b) 


Arguing  as  we  did  for  the  autoregressive  case,  no  a-priori 
knowledge  about  the  initial  values  of  the  moving  average 
coefficients,  implies  that. 
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[P(0|-1)] 


-1 


0. 


(2.35) 


And  even  if  after  we  processed  one  measurement  we  still 
admitted  no  knowledge  as  to  the  accuracy  of  the  previous 
estimate,  we  imply  that, 


[P(kjk-l)] 


-1 


0. 


(2.36) 


Substituting  these  implications  into  equations  (2.34)  our 
estimate  becomes. 


a(k+l|k)   =   [H(k)  r"""-  H'^(k)]"-'-  H(k)  r""^  z(k) 


(2.37) 


the  by  now  familiar  weighted  least  squares  estimate  [Ref .  28, 
29,30].   If  once  again  we  allow  R   to  be  unity,  let  k  tend 
toward  infinity  and  take  the  expectation  of  equation  (2.37), 
one  arrives  at  the  discrete  form  of  the  Wiener-Hopf  equation 
[Ref.  11] ,  namely 


a(k+l|k)   =   R   ~"^(k)  r   (k) 
—     '        uu       uz 


(2.38) 


where  it  can  easily  be  shown  that. 


R   (k; 
uu 


R   (0)    R   (-1) 
uu       uu 


R   (-1) 
uu 


R   (-m) 
uu 


R   (0) 
uu 


R   (-m) 
uu 

R   (1-m) 
uu 


R   (0) 
uu 


(2.39) 


and. 
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R   (k) 
uy 


R   (0) 
uy 

R   (-1: 
uy 


(2.40) 


R   (-m) 
uy 


Equation  (2.38)  provided  another  point  of  departure  from 
which  Perry  [Ref .  6]  develops  the  MA  Lattice  modeling  algorithms 
Appendix  A  develops  the  Wiener-Hopf  equation  otherwise  known 
as  the  all-zero  model  from  yet  another  approach. 

The  all-zero  model  is  fundamental  to  Widrow  least 
mean  square  (LMS)  adaptive  filters  [Ref.  11]  and  linear  pre- 
diction theory  [Ref.  16]. 

3 .   Autoregressive-Moving  Average  Form 

Returning  to  the  alternate  form  (equations  2.20)  of 
the  Kalman  Filter  equations  (2.19)  the  development  of  the 
Auto regressive  Moving  Average  (ARMA)  Adaptive  Kalman  Iden- 
tifier follows.   The  estimate  for  the  a.  and  b.  coefficients 

1      1 

is, 


a(k+l|k) 


b(k+l|k) 


[P(k+l|k)  [P(klk-l)  ]"^] 


a(klk-l) 


b(k|k-l) 


+  [P(k+l|k)  H'^(k)  R"-^]y(k) 


(2.41a) 


P(k+l|k) 


[  (P(k|k-1))"^  +H'^(k)  R"^H(k)]"-^ 


(2.41b) 
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Progressing  in  the  same  manner  as  in  the  previous  two  cases, 
we  assume. 


[P(0|-1)  ] 


-1 


0 


(2.42a) 


[P(klk-l)  ] 


-1 


0  . 


(2.42b) 


The   estimate    then   becomes. 


a(k+l   k) 

b(k+l   k) 

[H'^(k)  R  -^  H(k)]"-^  H(k)  R"-'-y(k)  .     (2.43) 


-1 


Taking  equation  (2.43)  a  step  further  by  letting  R    be 
unity,  letting  k  approach  infinity  and  taking  the  expectation 
we  have , 


a(k+l|k) 
b(k+l Ik) 


1-1  r 


R   (k) 

nu 


-R   (k-1) 
uv 


-R   (k-1)   R   (k) 
uy        yy 


R   (k) 

uy 


R   (k) 


(2.44) 


Note  that  the  time  varying  measurement  vector,  H(k),  is  of 
the  form 


H(k)   =   [u(k)  u(k-l)  . 


.  u(n-m)  I  -y(k-l)  -y(k-2: 


.. -y(k-n) ] 
(2.45) 


Shown  in  detail,  equation  (2.44)  has  the  characteristic  form 
(equation  2.46,  next  page).   It  is  the  Toeplitz  and  symmetric 
nature  of  equation  (2.46)  that  is  exploited  by  the  Levinson 
[Ref.  34]  and  Lattice  [Ref.  35]  algoritlims. 
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The  Adaptive  Kalman  ARMA  modeling  technique  can  be 
best  visualized  by  a  block  diagram.   Referring  to  Figure  2.1, 
the  unknown  system  is  excited  by  a  white,  zero  mean,  gaussian 
noise  sequence  of  sample  values  from  a  random  process.   The 
input  and  its  associated  output  are  then  passed  through  M 
and  N  delays  respectively  in  serial  form.   The  parallel  inputs 
(M+1)  and  the  outputs  (N)  are  concatenated  to  foirm  the 
measurement  vector  H(k)  which  is  represented  by  equation 
(2.13).   The  inputs  N  and  M  are  selected  a-priori  as  the 
model  orders  for  the  Autoregressive  and  Moving  Average  pro- 
cesses one  desires  to  identify.   It  can  be  easily  seen  that 
the  remainder  of  the  figure  simply  implements  equation  (2.15a) 
The  outputs  of  the  Adaptive  Kalman  Identifier  are  an  estimate 
of  the  coefficients. 


b(klk) 


a(k|k) 
and  a  one  step  prediction,  y(k|k). 

C.   OBSERVATIONS 

In  this  section  it  has  been  shown  that  a  direct  connection 
can  be  established  between  the  Adaptive  Kalman  Identifier  and 
the  Yule-Walker  equations  associated  with  an  AR  process. 
Secondly,  the  steady  state  Adaptive  Kalman  Identifier  closely 
resembles  the  discrete  Weiner-Hopf  equation  associated  with 
the  MA  process  when  the  measurement  matrix,  H(k) ,  is  of  the 
form, 
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H(k)   =   [u(k)  u(k-l)  u(k-2)  ...  u(k-m)] 

And,  thirdly  under  the  same  assumptions  as  in  the  previous 
two  cases,  the  Adaptive  Kalman  ARMA  Identifier  is  similar  to 
the  form  Perry  [Ref.  6]  exploits  using  Lattice  modeling. 
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III.   COMPARISON  BETWEEN  THE  ADAPTIVE  MA  KALMAN 

IDENTIFIER  AND  THE  WIDROW  LMS  ADAPTIVE  FILTER 

A.   PRELIMINARIES 

It  is  instructive  to  investigate  the  similarities  between 
the  Adaptive  Kalman  Identifier  and  the  Widrow  LMS  Adaptive 
Filter  when  the  concepts  are  applied  to  system  identifica- 
tion.  Basically,  the  LMS  algorithm  implementation  of  system 
identification  considers  a  block  diagram  as  is  shown  in 
Figure  3.1.   The  output,  y(k),  of  the  Adaptive  filter  is  simply 
a  weighted  linear  combination  of  the  past  and  present  known 
inputs.   The  same  input  is  fed  into  both  the  unknown  system 
to  be  identified  and  the  adaptive  filter.   The  output  of  the 
unknown  system  is  designated  the  desired  response,  d(k),  from 
which  an  error  signal,  e (k) ,  is  derived.   The  error  signal, 
£  (k) ,  provides  the  criterion  through  which  the  weights,  w., 
are  adjusted  such  that  the  error,  e (k+1) ,  is  driven  toward 
its  minimum.   A  more  detailed  analysis  of  the  operation  of 
the  LMS  algorithm  can  be  found  in  Ref.  11. 

The  weights  w.  are  adjusted  from  time  step  to  time 
step  in  the  following  manner. 


w(k+l)   =   w(k)  -  2k   £:(k)  X(k)  (3.1) 

—  —        s      — 

where  we  define  the  following  quantities. 
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w(k: 


WQ(k) 


w^(k) 
m 


(3.2) 


the  weight  vector  at  time  step  k. 


X(k)   = 


u(k) 
u(k-l) 


u (k-m) 


(3.3) 


the  input  signal  vector  at  time  step  k, 


e(k)   =   d(k)  -  y(k) 


(3.4) 


y(k)   =  w'(k)  X(k) 


(3.5) 


the  error,  £ (k) ,  and  the  output,  y(k),  at  time  step  k,  and, 
k  ,  a  scalar  constant  controlling  the  rate  of  convergence  and 
the  stability  of  the  adaptive  filter. 


B.   THE  COMPARISON 

Substituting  (3.4)  and  (3.5)  into  (3.1)  gives. 


w(k+l)   =  w(k)  -2k   X(k)  [d(k)  -w  (k)X(k)]. 


(3.6) 
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Recalling  the  MA  form  of  the  Adaptive  Kalman  Identifier, 
equation  (2.33a),  rewritten  here  for  convenience  as 

a(k+l|k)   =   a(k|k-l)  +K(k)  [z(k)  -H(k)  a(klk-l)],      (3.7) 

one  can  make  the  following  associations: 

a(k+l|k)  <=^  w(k+l)  (3.8a) 

a(klk-l)  <=^  w(k)  (3.8b) 

2(k)      <=±>  d(k)  (3.8c) 

H'^(k)     (z=^    X(k)  (3.8d) 

K(k)      <=^  -2k   X(j)  (3. Be) 

Recall   that   for   the  MA  form  of    the  AKI ,    the   measurement 
vector  H(k)    represents    a   vector   of   past   and    present    inputs. 
Namely, 

H(k)       =       [u(k)    u(k-l)     ...    u(k-m)]     .  (3.9) 

Therefore    the   associations    (3 . 8a) - ( 3 . 8d)    are    straightforward. 
However,    it   is   not   so   clear   as    to  what   is  meant  by   the 
association    (3.8e).      Digressing   a  moment   to   present  an  equi- 
valent expression    [Ref.    27]    for    the   Kalman   gain,    K(k), 

K(k)       =      P(k|k)    H*^(k)    r"-^  (3.10) 

and  substituting  (3.10)  into  the  association  (3.8e)  and  enter- 
taining the  conjecture  that  the  quantities  are  equivalent 
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under  certain  conditions ,  we  have : 

P{k|k)  H'^(k)  r"-"-   =   -2k   X(k)  .  (3.11) 

'  s  — 


The  conditions  alluded  to  are  (1)  that  the  Adaptive  Kalman 
Identifier  is  in  steady  state  and  (2)  that  the  statistics  of 
the  input  forcing  function  are  stationary.   Denoting  the  steady 
state  error  covariance,  PCkjk),  as  P^,  equation  (3.11)  can 
be  solved  for  k  , 

k  I   =   -  ^r"-^  P   .  (3.12) 

s         2      =° 

Invoking  the  entire  Kalman  gain  history  in  its  more  popular 
form  of  equation  (2.33b)  and  equating  it  to  the  Widrow  gain 
(3.8e),  we  instead  arrive  at, 


.2k^X(k)   =   P(klk-l)  HT(k)   [H(k)  P(k|k-1)  H^(k)  ^R]-^ 

(3.13a) 

H'^(k)   =   X(k)  .  (3.13b) 


Solving  for  the  convergence  factor,  k  ,  we  obtain, 

-2k  X(k)x'^(k)   =   P(klk-l)  X(k)  [x'^  (k)  P  (k  |  k-1)  X  (k)  +  R]  ""'■x'^  (k) 
-2k  I  =  P(k|k-l)X(k)  [x'^(k)P(k  lk-1)  X(k)  +R]"-^x'^(k)  [X(k)x'^(k)  ]"-^ 
But  [x'^(k)P(klk-l)X(k)  +R]"''"  is  a  scalar.   Therefore, 


P(k|k-1) 

-2k  I   =   -,= (3.14) 

X  (k)  P(k  ik-l)  X(k)  +  R 
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Griffiths  [Ref.  22]  suggests  that  the  convergence  factor  k  , 
which  he  denotes  as  |j(n),  should  be  chosen  such  that, 

M(n)   =   % (3.15) 

L  a^(n) 

X 

where  0  <  a  <  1  is  a  normalized  adaptive  stepsize  parameter 

^2 

and  the  term  a  (n)  is  an  estimate  of  the  input  power  level 

which  may  be  computed  using  a  geometrically-weighted  average 
for  an  L  weight  adaptive  filter. 


a^(n)   =   (1  -  ^)  Sj(n-l)  +  2L  x^  (n)  .  (3.16) 

X  Li       X  J-i 


Comparing  equation  (3.14)  and  (3.15),  one  observes  that  if 
P  (k  k-1)  is  equal  to  the  identity  matrix  then  we  obtain  a 
term  proportional  to  the  input  power.   Further  if  measurement 
noise  is  considered  negligible,  R  -^  0.0,  then  3.14  and  3.15 
are  essentially  equivalent.   The  salient  point  to  make,  however, 
is  that  the  choice  of  the  convergence  constant,  kg,  in  the 
LMS  algorithm  provides  no  clue  as  to  how  to  deal  with  measure- 
ment noise  whereas  the  Adaptive  MA  Kalman  Identifier  explicitly 
incoirporates  measurement  error  into  its  algorithm. 

C.   OBSERVATIONS 

It  has  been  known  that  the  LMS  algorithm  was  suboptimal 
since  the  actual  gradient  of  equation  (3.1)  is  replaced  by 
the  approximation,  2£  (k)  x(k),  [Ref.  24].   However,  the  degree 
of  suboptimality  is  difficult  to  quantify.   In  previous 
linear  prediction  research  no  mention  was  made  regarding  the 
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role  of  measurement  noise.   Equation  (3.14)  gives  us  some 
insight  into  wiser  selections  of  the  rate  of  convergence  con- 
stant, k  ,  since  it  takes  into  account  the  effects  of  measure- 
s 

ment  noise.   Further,  it  appears  that  the  LiMS  adaptive  filter 
is  a  special  case  of  the  Adaptive  MA  Kalman  Identifier. 
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IV.   COMPARISON  BETWEEN  THE  ADAPTIVE  ARMA  KALMAN  IDENTIFIER 
AND  THE  ADAPTIVE  RECURSIVE  L^.S  FILTER 

This  chapter  endeavors  to  investigate  the  similarities  and 
differences  between  the  Adaptive  ARMA  Kalman  Identifier  and 
the  Adaptive  Recursive  LMS  Filter  [Ref.  18].   LMS  filters 
[Ref .  11]  have  enjoyed  much  popularity  in  the  recent  past  due 
to  their  ease  of  implementation,  simple  unimodal  algorithm, 
robustness  and  ability  to  "adapt"  to  the  unknown  statistics 
of  the  signal  environment.   Being  transversal  in  nature,  they 
have  a  finite  impulse  response  being  able  to  produce  only 
zeros  in  the  input/output  transfer  function.   One  may,  how- 
ever, decide  to  model  the  transfer  function, 

H^(z)   = ^  (4.1) 

1  -  .  9z 

using  a  transversal  filter,  only  to  realize  that  a  large 
number  of  delays  are  required  in  order  to  arrive  at  an 
adequate  approximation.   The  germane  point  which  is  being 
made  is  that  one  weighted  feedback  tap  can  realize  an  infinite 
string  of  feed-forward  coefficients.   Moreover,  it  is  very 
desirous  to  adjust  the  feedback  weights  adaptively  in  some 
optimal  fashion  to  the  statistics  of  the  signal  environment. 

A.   PRELIMINARIES 

The  approach  taken  by  Feintuch  [Ref.  18]  is  patterned 
after  the  analysis  first  presented  by  Widrow  [Ref.  11].   A 
summary  of  Feintuch' s  derivation  will  be  presented  here  with 
emphasis  on  its  application  to  ARMA  modeling. 
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Recall  the  general  ARMA  equation  (2.1b), 


m  n 

y(k)   =    I      a.u(k-i)  +  I      b.y(k-i)  (4.2a) 

i=0   ^         1=1   ^ 


which  can  be   rewritten   in  vector  notation, 

y(k)       =      a\   +   b\  (4.2b) 


where , 


T 
a        =       [a_    a,     ...    a^]  (4.3a) 

—  0      1  m 


b*^      =       [b,    b_    .  ..    b    ]  (4.3b) 

—  12  n 

T 

Y_     =      [y(k-l)y(k-2)     ...    y(k-n)]  (4.3c) 


T 
u     =       [u(k)    u(k-l)     ...    u(k-in)]  (4.3d) 


Given  that  (4.2)  is  the  assumed  mathematical  description  of 
the  unknown  system  where  u  and  y,  the  input  and  output  data 
sequences,  are  known,  then  by  solving  (4.2)  for  a  and  b, 
identification  of  the  unknown  system  can  be  made.   The  LMS 
algorithms,  in  general,  employ  a  "desired"  signal  d(k)  with 
which  to  "train"  the  adaptive  filter.   If  the  desired  signal 
was  assumed  to  be  the  response  of  some  unknown  system  to  a 
known  input  signal,  then  the  algorithm  presented  by  Feintuch 
[Ref.  18]  can  be  used  as  a  means  by  which  to  identify  the 
system  parameters. 

The  first  step  in  the  Recursive  LMS  derivation  is  to  form 
the  error  between  the  desired  signal,  d(k) ,  and  the  output 
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of  the  filter,  y (k) , 

e(k)   =   d(k)  -  y(k)   =   d(k)  -  a'^u  -  b\  (4.4) 


Forming  the  square  of  (4.4)  and  taking  the  expectation, 

we  have  the  mean  square  error  representation  of  the  filter, 

E{e^(k)}   =   E{d^(k)}  +  a"^  R  a  +  b*^  R  b  ... 

—    uu—    —    yy— 

-  2a'^  R,   -  2b'^  R.   +  2a'^  R   b  (4.5) 

—   du    —   dy    —   uy— 

where , 

R   (k)   =   E{u(k)  u'^(k)}  (4.6a) 

R   (k-1)   =   E{y(k-1)  y'^(k-l)}  (4.6b) 


R,  (k)   =   E{d(k)  u(k)}  (4.6c) 

du 


R   (k)   =   E{d(k)  y(k) }  (4.6d) 

R   (k-1)   =   E{u(k-1)  y(k-l)}  .  (4.6e) 

J; 

It  is  assumed  that  the  statistics  (4.6)  are  constants  allowing 
the  gradient  of  (4.5)  to  be  taken  with  respect  to  the  a.  and 
b.  .   This  assumption  does  not  stretch  the  theory  since  in 
practice  one  uses  an  input  test  signal  with  stationary  char- 
acteristics and  if  the  unknown  system's  output  process  is 
not  stationary,  then  no  identification  of  the  system  parameters 
can  be  made.   The  respective  gradients  are. 
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|-[E{e2(k)}] 

a  di 


|^[E{e^(k)}] 


2R   (k)a  -  2R,  (k)  +  2R   (k)b 
uu    —    du       uy    — 


2R   (k)b-2R,  (k)  +  2R    (k)  a 
yy   —    dy        uy   — 


(4.7a) 


(4.7b) 


Since  the  second  order  statistics  are  assumed  to  be  known, 
equations  (4.7)  can  be  solved  for  a  and  b  by  setting  the 
gradients  equal  to  zero.   In  matrix  partitioned  form  we  have. 


R   (k)  !  R   (k) 

uu  I   uy 

R*^   (k)  !  R   (k) 

uy  I   yy   . 


R,  (k) 
du 


^dy(^) 


(4.8) 


At  this  point  we  digress  momentarily  to  mention  that  Johnson 
and  Larimore  [Ref.  19]  and  Widrow  and  McCool  [Ref.  20]  agree 
with  Feintuch's  derivation.   The  following  steps,  however, 
are  controversial  [Ref.  19,20]. 

Feintuch  continues  and  states  that  in  general,  the  statis- 
tics involved  in  (4.8)  are  not  available  a-priori .   One  method 
of  estimating  the  statistics  is  to  make  the  filter  adaptive 
in  an  iterative  fashion  using  the  method  of  steepest  descent. 
The  method  of  steepest  descent  employs  an  algorithm  of  the 
form, 


a(k+l)   =   a(k)  +k  |-[E{e^(k)}] 
—  —       a  oa 


(4.9a) 


b(k+l)   =  b(k)  +  k^   |^[E{e^(k)}] 


(4.9b) 


The  gradient  involved  in  (4.9)  can  be  approximated  using  the 
techniques  outlined  in  [Ref.  11].   The  iterative  LMS  algorithm 
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for  estimating  the  coefficients  a  and  b  is, 

a(k+l)   =   a(k)  -  2k  e (k)  u(k)  (4.10a) 

"~"  a 

b(k+l)   =   b(k)  -  2k^e(k)  y (k) .  (4.10b) 

Figure  4.1  is  a  block  diagram  implementation  of  the  Adaptive 
Recursive  LMS  Filter  algorithm  as  applied  to  system  identifi- 
cation.  The  input  signal,  u(k),  to  both  the  unknown  system 
and  the  adaptive  filter  is  a  zero  mean,  white  gaussian  noise 
sequence  of  samples  from  a  random  process.   The  output  of 
the  unknown  system  is  d(k)  which  is  used  as  the  training 
signal  for  the  filter.   The  output  of  the  adaptive  filter, 
y(k),  is  compared  with  the  desired  signal,  d(k),  from  which 
an  error  signal  is  derived.   The  error  signal,  e (k) ,  is  then 
used  in  the  LMS  algorithm,  equation  (4.10),  to  iteratively 
adjust  the  a  and  b  coefficients.   Theoretically,  when  the 
coefficients  of  the  zeros-producing  adaptive  filter  and  the 
poles-producing  adaptive  filter  correspond  exactly  with  those 
of  the  unknown  system,  the  error  sequence  should  be  zero  and 
the  unknown  system  is  modeled  by  the  Adaptive  Recursive  LMS 
Filter.   The  convergence  constants,  k   and  k,  ,  are  chosen 
a-priori  as  are  the  initial  values,  a(0)  and  b(0) ,  for  the 
coefficients . 

B.   THE  COMPARISON 

One  can  begin  the  comparison  between  the  Adaptive  Kalman 
Identifier  and  the  Adaptive  Recursive  LMS  filter  by  noting 
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the  similarities  between  the  equations  (2.43)  and  (4.8), 
repeated  here  for  convenience , 


a(k+l|k) 


b(k+l|k) 


R   (k) 
uu 


-R    (k-i: 
yu 


-R   (k-1) 
uy 


R  ^  (k) 

yy 


-1 


R   (k) 
uy 


R   (k) 

yy 


(2.43a) 


=   ^1  ^1 


(2.43b) 


R   (k) 
uu 


I  R   (k-1) 
I   uy 


R^   (k-1)  I  R   (k) 

uy    '  yy 


R,  (k) 
du 


R,  (k) 
dy 


(4.8a) 


R   (k) 
uu 


I  R   (k-i: 
I   uy 


1    -1 


R^   (k-1)  I  R   (k) 
uy       yy 


R,  (k) 
du 


W^ 


(4.8b) 


^2  S 


(4.8c) 


The  upper  left  partitions  of  B^  and  B_  are  obviously  identi- 
cal, in  that  the  AKI  and  the  Adaptive  Recursive  LMS  both 
implement  instantaneous  estimates  of  the  input  covariance 
functions,  R   (k) .   The  lower  right  partitions  of  B,  and  B^ 

uu  ■::)     f  12 

are  similar  with  one  subtle  difference.   Both  are  instantaneous 
covariance  functions;  however,  B,  employs  the  covariance, 
R   (k) ,  of  the  actual  data  whereas  B„  uses,  in  fact,  previ- 
ous filter  estimates  to  compute  the  output  covariance,  R   (k) , 
more  properly  denoted  by  R^'^  (k)  .   The  negative  sign  of  the 
upper  right  and  lower  left  partitions  in  B, ,  in  comparis 


on 


with  Bp,  is  a  result  of  the  initial  definition  for  the  recursive 
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weights,  b. .   That  is,  the  transfer  function  for  the  general 
ARMA  equation  can  be  written  as. 


«arma'^' 


m 


+     y     a.  z 

i=i   " 


-1 


-1 


m 

1  -    y    b. z 
i=i    ^ 


(4.9a) 


or, 


^'arma^^^  = 


m 


a-,  +   y   a.  z 
0    .  ^,   1 
1  =  1 


-1 


n 


(4.9b) 


1  +    y   b.  z 
i=i    ^ 


-1 


The  second  and  more  important  observation  is  that  cross- 
correlations  employed  by  equation  (4.8b)  use  the  instantane- 
ous values  for  the  past  estimates  of  the  outputs  while  (2.43a; 
uses  the  actual  past  values  for  y(k) .   Similarly,  comparing 
C,  with  C„,  one  finds  that  the  lower  partition  of  C„  uses  the 
cross-correlation  between  the  desired  signal  and  the  past 
estimates  of  the  adaptive  recursive  LMS  filter. 

Equation  (4.8a)  can  convey  the  above  information  more 

clearly  if  it  is  instead  written  as, 

-1 


R   (k) 
uu 


R   -(k-1) 
uy 


R  -(k-1) 
uy 


R-- (k) 

yy 


^■du^^: 


W' 


(4.10) 


The  convergence  factors,  k   and  k,  ,  can  be  compared  to 
the  steady  state  Adaptive  Kalman  Identifier  gains  by  making 
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the  identical  analogies  (3.8)  as  was  done  in  the  previous 
chapter.   One  also  observes  that  since  the  output  data  that 
is  processed  by  the  Adaptive  Recursive  LMS  filter  are  in 
actuality  filter  estimates,  a  more  correct  version  of  Feintuch's 
original  algorithm  (4.10)  is  written  as. 


a(k+l)   =   a(k)  -  2k  e (k)  u(k) 
—  —        a      — 


(4.11a) 


b(k+l)   =   b(k)  -  2k   e(k)  ^(k)  . 


(4.11b) 


In   agumented   form,    equation    (4.11)    is. 


a(k+l) 

a(k) 

2k  u(k) 
a   — 

b(k+l) 

b(k) 

_2kj^  y(k) 

a(k+l) 

a(k) 

2k  U(k) 
a   — 

b(k+l) 

■                                                     Mi 

b(k) 

2k^   y(k) 

[d(k)   -a(k)   u(k)  -   b(k)    y(k)  ] 


d(k)   -  [u(k)       y(k)  ] 


a(k) 

_b(k)^ 
(4.12) 


Recalling  equation    (2.15a)    in    a   slightly  different   form  we 
have , 


a(k+l  |k) 
b(k+l  Ik) 


a(k|k-l) 
b(k Ik-l) 


+   K(k) 


z(k)    -H(k) 


a(klk-l) 
b(k  ik-l) 


(4.13) 


It  is  a  simple  matter  to  explore  the  similarities  between 
(4.12)  and  (4.13)  and  make  the  following  associations. 
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AKI 


Adaptive  Recursive  LMS 


a(k+l|k) 


b(k+llk) 


a(k|k-l) 


b(klk-l) 


z(k) 


<=> 


<^ 


a(k+l) 
b(k+l) 

"a(k)" 
b(k) 


(=>   d(k) 


H(k)   =   [u(k)  ...     (=l)    [u(k)   y(k)] 
u  (k-m)   y(k-l) 
. . .  y (k-n) ] 


(4.14a) 


(4.14b) 


(4.14c) 


(4.14d) 


K(k) 


<=> 


2k   u(k) 
a 


}\  ■V(k)_ 


(4.14e) 


The  associations  (4.14d)  and  (4.14e)  are  the  major  differences 
between  the  AKI  and  the  Adaptive  Recursive  LMS.   It  has  been 
shown  that  the  mean  square  error  surface  employing  not  only 
estimated  feed  forward  coefficients  but  also  estimated  feed- 
back (recursive)  coefficients  is  in  general  multimodal  [Ref. 
36]  .   Hence,  the  Adaptive  Recursive  LMS  algorithm  does  not 
minimize  the  mean  square  error  [Ref.  19,20],  and  in  general 
the  gradient  algorithm  in  this  case  does  not  seek  the  global 
minimum. 
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C.   OBSERVATIONS 

The  complex  nature  of  estimating  recursive  coefficients 
is  yet  a  formidable  problem  [Ref.  19,20].   It  seems  that 
Feintuch's  algorithm  is  successful  due  to  the  a-priori  knowledge 
of  the  minimal  order  generating  the  desired  signal.   This 
knowledge  is  used  in  setting  the  order  of  the  adaptive 
recursive  filter. 
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V.   SYNTHETIC  DATA  GENERATING  PLANTS 

The  Adaptive  ARMA  Kalman  Identifier  was  tested  using 
computer  derived  data  from  several  models.   The  input/output 
data  was  collected  by  driving  a  known  plant  with  a  zero  mean, 
unit  variance,  white  gaussian  noise  sequence  of  samples. 
The  AKI  algorithm  was  repeatedly  tested  using  data  from 
progressively  more  complex  models. 

A.   AUTOREGRESSIVE-MOVING  AVERAGE  DATA 

Autoregressive-Moving  Average  (ARMA)  data  was  generated 
using  an  equation  of  the  form, 

y(k)   =   b,y(k-l)  +  b^y(k-2)  +  ...  +  b  y(k-n) 
"'  1  2  n 

+  a_u(k)  +  ...  +  a  u(k-m)  (5.1) 

U  m 

Taking    the    Z- transform  of    (5.1)    we   have, 

Y(z)       =      b,z~-^Y(z)    +   b^z~^Y(z)    +    ...    +   b    z~^Y(z) 
1  ^  n 

+   a^U(z)    +    ...    +   a    z""\j(z)  (5.2) 

0  m 

where  Y(z) (U(z))  represents  the  Z-transform  of  the  output 
(input)  and  z   Y(z)(z   U(z))  represents  the  Z-transform 
of  the  output  (input)  delayed  n.  time  steps.   Equation  (5.1) 
has  already  been  defined  as  the  general  form  of  the  ARMA 
process  provided  that  the  sequence  u(k-i)  i  =  0,1,...  comes 
from  a  gaussian  random  process  [Ref.  25],   The  ratio. 
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^    -1  ^        -m 
H<z,   =  ^4|I  =  ^^^1.1 ^B_ 

.  1  -  D,  Z    -  .  .  .  D  Z 

1  n 

can  then  be  readily  recognized  as  the  general  transfer  func- 
tion [Ref.  37]  of  a  digital  plant.   That  is,  the  u(k)  is  the 
observed  input  and  the  y(k)  is  the  perfectly  measured  output. 
Consider  that  the  measurements  of  the  output  include  some  ran- 
dom error,  v(k),  due  to  the  inaccuracy  of  the  measuring  instru- 
ments, then  the  resulting  situation  is  as  depicted  in  Figure 
5.1.   A  practical,  judicious  and  reasonable  description  of 
the  measurement  error  is  that  this  be  a  stochastic  process 
whose  distribution  is  zero  mean  gaussian.   The  noisy  measure- 
ment, z(k),  is  then  the  output,  y(k),  plus  the  measurement 
error,  v  (k)  . 

There  are  three  cases  of  interest:   Case  1:   b.  =  0  for 

1 

all  i;  Case  2:   a.  =  0,  i  =  1,2, ...,m;  Case  3:   a.  7^  0 , 

b.  7^  0  for  all  i. 

1 

1 .  Moving  Average  Data  (Case  1) 

Case  one  is  readily  recognized  as  the  all  zero  plant 
which  produces  moving  average  data.   A  simple  second  order 
plant  was  defined  where  a.  =  1.0,  a,  =  2.0  and  a^  =  3.0.   The 
gaussian  noise  sequences,  u(k)  and  v(k),  were  obtained  using 
the  general  purpose  IMSL  subroutine,  GGNML,  provided  by  the 
computer  center  at  the  Naval  Postgraduate  School. 

2.  Autoregressive  Data  (Case  2) 

Case  two  results  in  the  all  pole  plant  producing  auto- 
regressive data.   The  coefficients  were  selected  such  that  the 
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plant  of  Figure  5.1  produced  stable  data.   In  precise  terms, 
the  poles  of  the  function. 


H(z) 


1+24  ^   +24  ^   -^24^ 


(5.4) 


were  located  completely  within  the  unit  circle  in  the  Z-plane 
ensuring  a  stable  plant  [Ref.  37,38].   The  noise  sequences 
were  obtained  as  before. 

v(k) 


u(k) 


y(k) 


Figure  5.1.   ARMA  Digital  Plant 


3 .   Autoregressive-Moving  Average  Data  (Case  3) 

Case  three  follows  from  a  logical  combination  of  the 
two  previous  cases.   That  is,  the  transfer  function,  H(z),  now 
has  both  zeros  and  poles .   In  order  to  compare  with  previous 
work  done  in  the  area  of  system  identification  [Ref.  6],  one 
of  Perry's  models  [Ref.  6]  was  used.   Namely,  the  transfer 
function  H(z)  used  for  this  case  was. 


H(z)   = 


1  +  1.4z"-^  +  .98z"^ 


1  -  1.14Z"-'-  +  1.454 9 z"^  -  .88490z"^  +  .40745z"'^ 


(5.5) 
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The  input,  u(k),  was  obtained  as  in  the  previous  two  cases 
using  the  IMSL  subroutine  GGNML  set  to  unit  variance  and  zero 
mean.   The  output  measurement  noise  sequence,  however,  was 
set  to  a  variance  of  .0001.   The  basic  software,  program  ARMA, 
which  generated  the  different  data  is  included  as  Appendix  B. 
Basically,  the  coefficients  Al-AlO  change  the  character  of  the 
general  ARMA  equation  (5.1)  which  can  be  selected  to  produce 
the  desired  plant  data.   The  input/output  data  is  then  written 
onto  a  disk  file  for  subsequent  analysis. 

In  order  to  compare  some  of  the  findings  in  this  the- 
sis with  the  results  previously  obtained  by  Feintuch  [Ref.  18], 
the  transfer  function, 

zjts                         0.05  -  0.40z~''"  ,^  ^> 

H(z)   =   ^ J  (5.6) 

1.0  -1.1314z"-^  +0.25z" 

was  also  used  as  a  source  of  synthetic  data.   Feintuch  used 
(5.6)  as  a  source  of  data  with  which  to  analyze  the  operation 
of  his  Adaptive  Recursive  LMS  Filter. 

B.   PHASE  LOCKED  LOOP  DATA 

The  second  class  of  data  used  to  exercise  the  Adaptive 
ARMA  Kalman  Identifier  was  derived  from  a  computer  simulated 
phase  locked  loop  (PLL)  developed  by  Romeo  [Ref.  7].   The 
basic  PLL  algorithm,  however,  implemented  a  forward  Euler 
integration  scheme. 

Briefly,  in  block  diagram  form,  the  phase  tracking 
characteristics  of  the  PLL  in  the  frequency  domain  can  be 
depicted  as  in  Figure  5.2. 
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Y(s) 


Figure  5.2.   PLL  Block  Diagram 


Where : 


A   is  the  input  voltage 
K   is  the  loop  gain 
F(s)   is  the  filter  transfer  function 


For  most  applications,  it  is  assumed  that  cj)(t)  <<  77/2  allowing 
one  to  make  the  approximation,  sin  ex  z   a   .       This  is  the 
linear  mode  of  operation.   Romeo  [Ref.  7]  chose  the  filter 
characteristics  F(s)  =  1  +  K/S  and  the  same  characteristics 
were  used  in  this  thesis  in  order  to  compare  results.   The 
parameters  of  the  overall  system  were  adjusted  to  obtain 
a  damping  coefficient,  ;,  of  0.3  and  a  natural  frequency,  w 
of  3.33  rad/sec.   This  resulted  in  a  step  response  overshoot 
of  about  46.7%  at  t  ~  .75  sec.   Solving  for  the  time  domain 
step  response  analytically  in  the  linear  region  one  obtains. 


y(t)   =   1  +  1.048e"^  sin (3 . 178t  -  1 . 266) 


(5.7) 


The  final  block  diagram  of  the  PLL  system  is  shown  in  Figure  5,3 


U(S) 


kI 

Sin(-) 

2.0 

K  = 
5.55 

1 
s 

Y(s) 


Figure  5.3.   Final  PLL  Block  Diagram 

Normalized  step  responses  for  several  step  magnitudes 
using  Digital  Simulation  Language  (DSL)  were  obtained.   Figure 
5.4  shows  the  step  responses  of  several  step  inputs  in  the 
range  0  <  u(t)  <  30.   It  is  obvious  that  the  PLL  is  in  its 
approximate  linear  region  and  does  not  exhibit  any  discernible 
distortion.   Figure  5.5  is  a  repeat  of  the  above  test;  how- 
ever, the  range  of  step  inputs  are  30  <_  u(t)  <_   170.   It  is 
apparent  from  Figure  5.5  that  the  sin  (•)  nonlinearity  begins 
to  have  some  effect  on  the  operation  of  the  PLL  in  that  the 
amplitude  of  the  responses  are  reduced  and  delayed. 

The  DSL  PLL  system,  however,  was  not  used  as  a  source  of 
synthetic  data  because  of  the  problem  of  nonstationary  statis- 
tics of  the  input  signal  using  the  random  sequence  generators 
available  under  DSL  already  discussed  in  Ref.  7.   These 
simulations  were  nevertheless  used  as  a  basis  for  comparison 
of  the  operation  of  the  discrete  PLL  developed  in  Reference  7 
but  modified  to  implement  a  forward  Euler  integration  scheme 
which  is  used  in  this  thesis. 
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2ND  ORDER  PLL  --  TRflPEZQIDRL  INTEGRflilQN 
NORMRirZED  STEP  RESPONSE  . 00 1 /. 0 1 /,  1 /  1  / 1 0/20/30 


l.QC 


2.Q0 


3.00 


4.00 


5.00 


6.00 


XSCflLE=  1.00 
YSCflLE=  0.20 


UNITS/INCH 
UNITS/INCH 


RUN      NO.    1 
PLOT    NO.    1 


Figure    5.4.      PLL   Step    Response,    Linear    Region 
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2ND  ORDER  PLL  --  TRflPEZaiDRL  INTEGRRTIQN 
ORMRLIZED  STEP  RESPONSE  1 /30/60/9G/ I  20/ 1 50/ 1  70 


1.00 


4.00 


8,00 


9.00 


XSCflLE=  l.QO 
TSCRLE=  Q.20 


UNITS/INCH 
UNITS/INCH 


RUN      NO.    1 
PLOT    NO.    1 


Figure    5.5.      PLL   Step    Response,    Non-Linear   Region 
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The  forward  Euler  integrator  takes  the  form. 


i<=^ 


Tz 


-1 


1  -  z 


-1 


T  =  .01 


(5.8) 


The  block  diagram  of  Figure  5.3  changes  character  to  appear 
as  in  Figure  5.6. 


U(z)  ^ 


^ 

'  + 

> 

1 

5in(-) 

2.0 

K  = 
555 

-^ 

Tz-1 

Tz-^ 

Y(z) 


Figure  5.6.   Discrete  PLL  Block  Diagram 

Applying  Mason's  gain  rule  to  the  block  diagram  of  Figure  5.6 
in  the  linear  region  we  arrive  at, 


Y(z) 
U(z) 


=  H(z; 


0.02z"-^  -  0.01889Qz"^ 
1  -  1.980z"-^  +  .981110z"^ 


(5.9) 


the  linearized  PLL  transfer  function.   The  discrete  PLL  sys- 
tem has  zeros  at  (0.0,. 9445)  and  a  complex  conjugate  pole 
pair  at  (0.990  ±  j. 03178)  . 

The  discrete  PLL  system  was  implemented  using  a  Fortran 
program  on  the  IBM  370  in  double  precision.   The  source  code 
for  the  discrete  PLL  is  included  as  Appendix  C.   The  discrete 
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PLL  was  tested  using  various  magnitude  step  inputs.   Results 
were  similar  to  those  from  the  DSL  simulation.   That  is,  the 
discrete  PLL  exhibited  the  same  approximate  overshoot,  rise 
time,  natural  frequency  and  nearly  identical  time  of  max  over- 
shoot in  the  linear  region.   To  implement  the  PLL  using  forward 
Euler  integration,  the  block  diagram  of  Figure  5.6  was  recon- 
figured to  appear  as  Figure  5.7.   Note  that  Figure  5.7  does 
not  employ  any  delay  free  loops,  hence  it  is  easily  programmable 


u(k-l) 


e(k-l) 


Sin(-) 


v(k-l) 


,02 


1-.9445Z 


-1 


l-2z  Kz  ^ 


y(k) 


Figure  5.7.   Programmable  PLL  Implementing 
Forward  Euler  Integration 


Figure  5.8  summarizes  the  results  of  the  discrete  PLL  simula- 
tion for  normalized  step  inputs  in  the  range  0  <  u(t)  <_   30. 
Similarly,  the  discrete  PLL  was  tested  in  the  nonlinear  region 
using  normalized  step  inputs  in  the  range  30  <_   u(t)  <_  170. 
The  discrete  PLL  displayed  the  same  characteristics  of  de- 
creased output  amplitude  and  delay  which  was  first  observed 
in  the  DSL  simulation.   Figure  5.9  summarizes  the  test  results 
of  the  discrete  PLL  in  the  nonlinear  region. 
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Figure  5.8.   Discrete  PLL  Step  Response,  Linear  Region 
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Figure  5.9.   Discrete  PLL  Step  Response,  Nonlinear  Region 
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It  was  felt  that  the  PLL  was  adequately  modeled  for  use 
as  a  source  of  data  to  be  subsequently  analyzed  by  the  Adap- 
tive Kalman  Identifier. 
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VI.   MODEL  IDENTIFICATION  SOFTWARE 

Three  basic  Fortran  IV  programs  for  use  on  the  IBM  370 
were  developed: 

(1)  ADAPTSN:   Adaptive  Kalman  Identifier 

(2)  LMS :   Adaptive  Least  Mean  Square  Filter 

(3)  LMSR;   Adaptive  Recursive  Least  Mean  Square  Filter 

All  three  programs  used  double  precision  arithmetic  to  minimize 
the  effects  of  truncation  error,  limit  cycles  and  roundoff 
error.   In  all  cases  the  orders  of  the  autoregressive  and 
moving  average  processes  were  read  in  unformatted  form  from 
a  disk  file  along  with  the  input/output  data  of  the  unknown 
system  to  be  analyzed.   An  attempt  was  made  to  minimize  memory 
storage,  but  not  at  the  expense  of  program  flow  and  clarity. 

A.   ADAPTIVE  KALMAN  IDENTIFIER 

Program  ADAPTSN  is  a  versatile  Fortran  IV  software  program 
which  implements  equations  (2.15)  to  identify  the  coefficients 
of  an  MA,  AR,  or  ARMA  process.   Additionally,  at  each  iteration 
the  poles  and  zeros  of  the  evolving  transfer  function  are  com- 
puted.  ADAPTSN  can  also  be  used  to  perform  regression  analysis 
on  non-linear  terms  as  discussed  in  Chapter  VII.   The  versa- 
tility of  the  AKI  lies  in  the  various  options  which  are  avail- 
able to  the  user.   By  properly  selecting  the  flags  N,  M,  and 
NL,  the  user  can  perform  regression  analysis  on  data  whose 
combined  order  (N+M+1)  is  less  than  or  equal  to  20.   Table  6.1 
lists  the  ascribed  meanings  of  the  various  flags. 
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Table  6.1 
AKI  FLAG  MEANINGS 

N:   the  order  of  the  autoregressive  process 

M:   one  greater  than  the  order  of  the  moving  average 
process 

NL  =  0 :   AKI  is  used  to  identify  the  coefficients  associated 
with  the  general  ARMA  equation  (2.1b)  (linear 
regression) 

NL  =  1 :   AKI  is  used  to  identify  the  coefficients  of  the 
linear  general  ARMA  equation  and  the  weighting 
coefficient  associated  with  u-^(k-i) 

NL  =  2:   AKI  is  used  as  in  NL  =  0,1  and  additionally  iden- 
tifies the  weighting  coefficient  associated  with 
y3(k-i) 

NL  =  3 :   AKI  is  used  as  in  NL  =  0,1,2  and  additionally  iden- 
tifies the  weighting  coefficient  associated  with 
u2(k-i)y(k-i) 

NL  =  4 :   AKI  is  used  as  in  NL  =  0,1,2,3,  and  additionally 

identifies  the  weighting  coefficient  associated  with 
y2  (i-i) u (k-i) 

NL  =  5 :   AKI  is  implemented  to  analyze  time  series  and 
identify  the  ARMA  coefficients  associated  with 
the  Box-Jenkins  model  [Ref.  26]. 


The  option  given  by  NL  =  5  was  not  extensively  tested  and  as 
such  only  limited  results  are  available.   Table  6.2  outlines 
the  allowable  flag  combinations  which  will  produce  a  valid 
analysis  of  the  data. 

For  purposes  of  this  thesis,  the  non-linear  (NL)  options 
were  configured  to  implement  the  expansion  terms  associated 
with  the  Taylor  series  expansion  of  the  sine  function.   This, 
however,  can  be  changed  at  the  user's  discretion  to  implement 
any  other  series  expansion  by  inserting  the  appropriate  Fortran 
statements  at  the  proper  location  in  the  AKI  program. 
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Table  6.2 
VALID  FLAG  COMBINATIONS 

OPTION  (NL)  MAX  COMBINED  ORDER 


0  N  +  M  <_  20 

1  N  +  2M  <_  20 

2  2N  +  2M  <_   20 

3  3N  +  2M  ^  20 

4  4N  +  2M  ^  20 

5  N  +   M  <  20 


N,M  7^  0 

N,M  7^  0 

N,M  7^  0 

N,M  7^  0 


The  overall  structure  of  the  AKI  program  uses  subroutine 
calls  to  compute  the  various  quantities  necessary  for  the 
eventual  computation  of  the  system  ARMA  coefficients.   These 
subroutines  are:   GAIN,  RECKON,  PRINT,  NEXT.   A  brief  des- 
cription of  the  function  each  performs  is  given  at  the 
beginning  of  each  subroutine.   Figure  6.1  describes  the 
general  program  flow  of  the  AKI .   The  source  code  for  program 
ADAPTSN  is  included  as  Appendix  D. 

B.   ADAPTIVE  LEAST  MEAN  SQUARE  FILTER 

Program  LMS  uses  the  equations  developed  in  Chapter  III 
to  compute  the  coefficients  (or  weights)  associated  with  a 
moving  average  process.   The  LMS  filter  program  is  capable  of 
computing  up  to  12  MA  coefficients;  or  equivalently ,  is  capable 
of  identifying  an  eleventh  order  moving  average  process .   The 
general  program  flow  of  the  LMS  filter  is  presented  in 
Figure  6.2  and  the  annotated  source  code  is  included  as  Appen- 
dix E. 
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(2)  data 
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Figure    6.1.      AKI   General   Program   Flow 
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Call  GAIN: 
Caipute  G(k) 


Eon.    (2.15b) 
Call  RECKON: 


conputB  a, 6 
Eqn    (2.15a) 


cnrprcs 

P(k+llk), 
P(klk) ,   Eq. (2. 


ISt) 
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CCrpute  Poles 
and  Zeros 


Call  NEXT 


Figure    6.1.        (CONTINUED) 
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Figure    6.2.      LMS   Filter   Program  Flow 
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C.   ADAPTIVE  RECURSIVE  LEAST  MEAN  SQUARE  FILTER 

Program  LMSR  realizes  the  algorithm  proposed  by  Feintuch 
and  recapped  in  Chapter  IV.   It  was  designed  to  handle  a 
combined  AR  and  iMA  order  of  11.   That  is,  the  program  can  be 
used  to  estimate  m  +  1  =  M  coefficients  associated  with  the 
MA  process  and  n  =  N  coefficients  associated  with  the  AR 
process  such  that  N  +  M  ^  12.   The  general  program  flow  for 
the  Adaptive  Recursive  LMS  filter  is  shown  in  Figure  6.3  and 
the  source  code  is  included  as  Appendix  F. 
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Figure  6.3.   LMSR  Filter  Program  Flow 
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VII.   NON-LINEAR  IDENTIFICATION 

The  Adaptive  Kalman  Identifier  can  be  modified  to  estimate 
the  weighting  coefficients  associated  with  higher  order  terms 
produced  by  some  non-linearity  within  the  unknown  system.   In 
general,  much  has  to  be  known  about  the  non-linearity  such 
that  the  functional  description  chosen  for  it  is  a  close 
approximation  to  the  effects  it  causes.   In  this  thesis  one 
possible  application  of  the  Adaptive  Kalman  Identifier  toward 
identifying  a  system  with  a  known  non-linearity  is  explored. 
No  attempt  has  been  made  to  present  an  extensive  treatment  of 
non-linear  analysis  techniques. 

The  approach  taken  has  been  previously  explored  by  Parker 
[Ref.  8].   A  special  case  of  the  generalized  non-linear  ARMA 
model  [Ref .  8] , 

oo  oo  oo 

y(k)    =        I    a(i,)u(k-i,)   +       I  [    a (i, ,i^) u (k-i^ ) u (k-i^)   +... 
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00  oo 

+        I  ...         5[   b(j^,  j^,  .  .  .  ,  j^)y(k-j^)y  (k-j^)     ...y(k-j^) 

-"l  -"n 

30  00                                                                                                                                                                                        00                                         OO                        OO 

+        I  j;   c(i    ,j    )u(k-i    )y(k-j    )   +  ...  +       I      -"      I           I 

i  .=0  j    =1        ^      ^               ^               ^                    i   =0           i   =0    j    =1 

]  -'I                                                                             1                  m        -"1 


...       I    c(i-|_  /  .  .  .i^/ j-j^.  .  .  j^)  u(k-i^)  .  .  .u(k-i^)  y  (k-jj^)  .  .  .y  (k-j    ) 

(7.1) 


^n=l 
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is  used  to  model  the  input/output  relationships,  non-linear 
transfer  function,  of  the  phase  locked  loop.   The  non-linear 
element,  the  sine  function,  of  the  PLL  can  be  replaced  by 
a  Taylor  power  series  expansion, 


sin(x)   =   X  -  YT  ^   ■*■  "ET  ^ 


(7.2) 


Other  expansions  can  be  used,  e.g.,  Legendre  polynomials, 
Volterra  series,  ...,  etc.;  however,  only  the  Taylor  expan- 
sion was  investigated  in  this  thesis.   Practical  implementa- 
tion of  (1.2)    suggests  truncation  at  the  third  order  term. 
The  sine  function  is  therefore  approximated  as 


sin (x) 


1   3 
X  -  -  X 


(7.3) 


Substituting  (7.3)  into  the  block  diagram  of  Figure  (5.6)  we 
have  Figure  7.1. 


z"\(z) 


z"-^E2(z) 


z  u(z)  +, 


^0- 


,-i<..^ 


.02 


1-.9445Z 


-1 


l+2z  ■^+z  ^ 


Y(z) 


Figure  7.1. 


PLL,  Third  Order  Taylor  Series 
Approximation 
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The  discrete  time  domain  equations  at  the  different  nodes  are 


e^(k-l)   =   u(k-l)  -  Y(k-l)  (7.4a) 

62(^-1)   =   e^(k-l)  -j^  ej(k-l)  (7.4b) 


y(k)   =   .02e  (k-1)  -  .01889e  (k-2)  +  2.0y(k-l) 

-  y(k-2)  (7.4c) 

Manipulating    equations     (7 . 4a) - (7 . 4c)    gives 

y(k)       =      .02u(k-l)    -  .  01889u  (k-2)   +1.98y(k-l)    -  .  98111y  (k-2) 

-  .003333u^ (k-1)    +    .003148333u^(k-2) 

+  .003333y^(k-l)  -  . 003148333y^ (k-2)  (7.5) 

+  .Olu^ (k-l)y(k-l)  -  .09445u^(k-2)y(k-2) 

-  .01u(k-l)y^ (k-1)  +  .09445u(k-2)y^ (k-2) . 

Equation  (7.5)  indicates  which  non-linear  terms  of 
equation  (7.1)  should  be  retained.   Therefore,  the  third 
order  non-linear  approximation  model  should  contain  the 
following  terms: 

u(k-l,  u(k-2),  y(k-l),  y(k-2)      linear  terms         (7.6a) 

3        3 
u  (k-1),  u  (k-2)  input  cubic  terms    (7.6b) 

3        3 
y  (k-1),  y  (k-2)  output  cubic  terms   (7.6c) 
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u^(k-l)y(k-l)  ,  u^(k-2y(k-2)   ]  ^v,  •  ^    ^ 

^    -^             ^  third  order  cross    ,.,  ^-,, 

\.  J      J.    J.                             (7 .  6d) 

_              9             1  product  terms 

y^(k-l)u(k-l) ,  y^(k-2)u(k-2)   J 

The  AKI  algorithm  is  primarily  modified  to  include  the 
hybrid  signals  (7.6)  in  the  structure  of  H(k),  the  measure- 
ment vector.   Since  now  the  measurement  vector,  H(k)  takes 
the  form, 

H(k)  =  [u(k)  ...  u(k-m),-y(k-l)  ...  -y(k-n)  ,u^(k)  ... 

...  u^(k-m),-y^(k-l)  ...  -y^(k-n)  ,-(u^  (k-l)y(k-l) )  . . . 

...-(u^(k-n)y(k-n)),(u(k-l)y^(k-l))  ...  (u  (k-n)  y^  (k-n) )  ] 

(7.7) 

the  AKI  algorithm  calculates  the  coefficients  associated  with 
the  special  case  of  the  generalized  non-linear  ARMA  model, 
equation  (7.1) . 

Romeo  using  a  similar  approach  [Ref.  7]  computes  a  least 
squares  curve  fit  for  the  third  order  truncation  of  the 
Taylor  series  for  the  sine  function  over  the  interval 
(0,7t/2).   Beginning  with  the  approximation, 

3 

sin(x)  z      ax  +  3x  (7.8) 

Romeo  finds  the  values  for  a  and  S  to  be  .865  and  -.095 
respectively.   Perfoinning  the  same  operations  as  those  used 
in  obtaining  equation  (7.5)  we  arrive  at. 
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y(k)  =  .02au(k-l)-.01889au(k-2)  +  (2.0  -  .02a)y(k-l) 

+  (.01889a -l)y(k-2)  -  .026u-^(k-l)  +  .018896u^(k-2) 

-  .0188937"^  (k-2)  +  .066u^  (k-l)y  (k-1)  +  .026y^(k-l) 

-  .0566706u^(k-2)y(k-2)  -  .066y^(k-l)u(k-l) 

+  .056670ey^(k-2)u(k-2)  (7.9) 


Table  7.1  tabulates  the  resulting  weighting  coefficients 
using  the  three  methods  just  described:   (1)  analytically 
calculated  values  using  the  truncated  Taylor  series,  (2)  least 
squares  estimate  of  the  third  order  sine  approximation,  as 
computed  by  Romeo  [Ref.  7]  and  (3)  the  coefficient  estimates 
as  computed  using  the  AKI  algorithm. 

The  tabulation  clearly  shows  that  the  AKI  outperforms 
the  first  two  methods  overall.   That  is,  the  linear  terms 

were  identified  without  question;  however,  the  AKI  failed 

3 
to  identify  the  coefficients  associated  with  the  y  (i-1) , 

3         2  2 

y  (k-2) ,  u  (k-2)  y(k-2)  and  y  (k-2)  u(k-2)  terms.   The  reason 

for  the  failure  is  not  known  and  was  not  investigated. 
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Table  7.1 


NON-LINEAR  WEIGHTING  COEFFICIENTS  CALCULATED 
USING  THREE  METHODS 


(1) 


TAYLOR  SERIES 
Sin  (x)   =  X  -  -jj-x 


y(k-l) 

y(k-2) 

u(k) 

u(k-l) 

u(k-2) 

u^k) 

u^(k-l) 

u^  (k-2) 

y^k-1) 

y-^(k-2) 

u  (k-l)y(k-l) 

u^{k-2)y(k-2) 

y^k-l)y(k-l) 

y^  (k-2)  u  (k-2) 


1.98 

-.98111 

0.0 

0.2 
-.01889 
0.0 
-.003333 

.003148333 

.003333 
-.003148333 

.01 
-.09445 
-.01 

.09445 


(2) 

LEAST  SQUARES  EST. 
sin  (x)  =  ox  +  Sx^ 
a  =  .865,  3  =  -.095 

1.982700 

-.98366015 

0.0 

.01730000 
-.01633985 
0.0 

.00190 
-.00179455 
-.00190 

.00179455 
-.005700 

.00538365 

.005700 
-.00538365 


(3) 

AKI  USING  ACTUAL 
sin(x)  data  base 


1.980 
-.9815 

.00003691 

.02003 
-.01892 

.00001471 
-.003278 

.003032 
-.05695 

.05799 

.01032 
-.009610 
-.01640 

.01913 
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Vin.  FINDINGS  AND  CONCLUSION 

The  performance  of  the  Adaptive  Kalman  Identifier  was 
compared  to  the  LMS  Adaptive  Filter  and  the  Adaptive  Recur- 
sive LMS  Filter  using  data  derived  from  the  models  discussed 
in  Chapter  V.   Results  for  the  case  where  the  order  of  the 
model  (that  is,  m  and  n  of  equation  (2.1b))  are  assumed 
known  are  presented  first  followed  by  the  analysis  of  the 
PLL  data.   The  findings  for  the  overmodeled  case  are  presented 
next.   The  graphs  show  typical  runs  and  do  not  represent 
ensemble  averages. 

A.   ORDER  OF  THE  UNKNOWN  SYSTEM  IS  KNOWN 
1.   AKI  vs.  LMS  Adaptive  Filter 

Using  the  MA  form  of  the  Adaptive  Kalman  Identifier, 
its  performance  can  be  compared  against  that  of  the  LMS 
adaptive  filter.   Synthetic  data  derived  from  a  plant  whose 
transfer  function  is, 

H(z)   =   1.0  +  2.0z"-'-  +  3.0z"^  (8.1) 

was  used.   However,  a  fair  comparison  necessitated  that  the 
LMS  Adaptive  filter  be  "tuned"  by  adjustment  of  the  conver- 
gence factor,  k  .   Several  values  for  k   in  the  range  [-.600, 
-.200]  were  used.   The  objective  of  tuning  the  LMS  filter  was 
to  achieve  a  fast  convergence  time  with  little  or  no  steady 
state  error  while  not  compromising  filter  stability.   The 
three  filter  weights  were  normalized  and  plotted  for  five 
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convergence  factor  values.   It  can  be  seen  from  Figures  8.1- 
8.3  that  convergence  is  essentially  reached  by  step  thirty- 
five,  (k  ~  35).   As  the  convergence  factor  is  further  in- 
creased it  can  be  noted  in  Figure  8.4  that  the  filter  weights 
become  more  noisy,  and  that  when  k   =  -.600,  filter  stability 
is  being  compromised  (Figure  8.5a). 

Using  the  same  data,  the  AKI  converges  to  the  iMA  coeffi- 
cients in  less  than  five  iterations.   Referring  to  Figure 
8.5b,  it  was  also  noted  that  as  the  measurement  noise  v(k) 
was  increased,  the  LMS  algorithm  yielded  more  noisy  estimates 
whereas  the  AKI  tended  to  compensate  for  measurement  noise. 
This  is  intuitively  reasonable  since  the  AKI  incorporates 
into  its  algorithm  the  effects  of  measurement  error  due  to 
noisy  sensors.   Comparing  Figure  8.6  with  Figure  8.7  the 
latter  figure  clearly  shows  the  faster  convergence  of  the  AKI 
to  the  plant  coefficient. 

It  was  shown  in  Chapter  III  that  the  AKI  gains  would 
approach  a  steady  state  value  and  indeed  they  did,  as  Figure 
8.8  shows.   However,  the  instantaneous  Kalman  gains  display 
a  more  erratic  pattern  as  shown  in  Figure  8.9.   The  informa- 
tion of  Figure  8.8  was  gleaned  from  Figure  8.9  by  computing 
the  average  of  the  individual  gains  at  each  iteration  using 
the  algorithm, 

E{G(k+l)}   =   E{G(k)}  +  kTTtG(k)  -E{G(k)}]  (8.2) 

whe  re , 


E{G (k+1) }   is  the  one  step  prediction  of  the 
average  value  for  the  gains. 
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Figure  8.1.   LiMS  Adaptive  Filter  Weights,  k^  =  -.200 
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H(z)  =1.0+2.0z~-^  +  3.02~'^ 


DISCRETE    TIME,    K 


Figure  8.2.   LMS  Adaptive  Filter  Weights,  k   =  -.300 
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H(z)   =  1.0+2.0z~-^  +  3.0z~^ 


DISCRETE    TIME.    K 


Figure  8.3.   LMS  Adaptive  Filter  Weights,  k   =  -.4  00 
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Figure  8.4.   LMS  Adaptive  Filter  Weights,  k   =  -.500 
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Figure  8.5a.   LMS  Adaptive  Filter  Weights,  k   =  -.600 
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Figure  8.7.   LiMS  Adaptive  Filter  Weights 
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Figure  8.8.   AKI  Averaged  Gains  for  Moving  Average  Model 
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E{G (k) }   is  the  present  average  value  for  the  gain 

G(k)    is  the  present  value  of  the  instantaneous 
gains 

T-^^       is  the  averaging  gain  . 

This  method  of  obtaining  a  "running"  average  is  generally 
well  known,  (see  for  example  Ref.  39).   The  example  presented 
was  not  the  only  one  used  but  a  good  representative  of  the 
operation  of  the  AKI  when  the  data  was  derived  from  a  moving 
average  process. 

2.   AKI  Applied  to  Autoregressive  Data 

Before  applying  the  AKI  to  identifying  the  coefficients 
of  a  complex  ARMA  plant,  it  was  first  tested  using  data 
derived  from  a  plant  whose  transfer  function  specified  an 
autoregressive  process.   The  transfer  function  used  was. 


H(z) 


,  ^  ^  26.0  -1  ^  9.0  -2  ^  1.0  -3 


^  (8.3) 


l-b,z    -bz    -b-,z 
12        3 


The  AKI  had  no  problem  converging  on  the  plant  coefficients 
(Figure  8.10)  producing  the  following  estimates  at  the  4  2 
iteration  (Table  8.1).   It  essentially  converged  on  the  plant 
coefficients  in  approximately  five  iterations.   The  gains 
computed  by  the  AKI  for  this  case  also  displayed  convergence 
to  a  steady  state  value.   The  gain  history  for  this  example 
is  depicted  in  Figure  8.11. 
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Figure  8.11.   AKI  Autoregressive  Averaged  Gains 
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Table  8.1 
ACTUAL  VS  AKI  ESTIMATES  FOR  COEFFICIENTS  OF  EQUATION  8.3 

ACTUAL  AKI  ESTIMATES 

dQ      1.0  0.9995 

b,     -1.08333  ...  -1.083 

h^  -0.3750  -0.3756 

b^     -0.04166  . . .  -0.04212 

3.   AKI  Applied  to  ARMA  Data 

The  next  logical  step  was  to  use  the  AKI  to  identify 
the  coefficients  of  a  general  ARMA  process.   One  of  Perry's 
models  [Ref.  6]  was  used  for  this  purpose.   Specifically,  the 
transfer  function  of  the  plant  was, 

„,  ,  1.0  +  1.4z"-^  +  .98z^ 

H(z)   =   zi r2 =3 =4 

1-1.142  ^+1.45492   -.884902   +.40745z 

^0  ^  ^1^    ""  ^2^ 
^  ^  ^  (8.4) 


1  -  b,2    -  b^2    -  b-,z    -  b.z 
12        3        4 


The  ARMA  plant  was  subjected  to  the  same  conditions  which 
Perry  [Ref.  6]  describes.   That  is,  unit  variance,  zero  mean, 
white  gaussian  noise  was  used  as  the  input.   The  reader  is 
reminded  that  the  output  signal  processed  by  the  AKI  was  cor- 
rupted by  measurement  noise  as  described  in  Chapter  V.   The 
output  data  used  in  Perry's  examples,  however,  reflects  noise- 
less measurements.   Table  8.2  tabulates  the  results  for  the 
coefficient  estimates  computed  by  the  AKI .   Even  though  the 
results  presented  in  Table  8.2  represent  the  coefficient 
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ACTUAL 

^0 

1.00000 

^1 

1.40000 

^2 

0.98000 

^3 

0.00000 

^4 

0.00000 

\ 

1.14000 

^2 

-1.45490 

'=3 

0.88490 

"^4 

-0.40745 

Table  8 . 2 

ACTUAL  VS  AKI  ESTIMATES  FOR  COEFFICIENTS  OF  EQUATION 
8.4  AT  THE  371ST  ITERATION 


AKI  ESTIMATES 

0.98760 

1.40100 

0.99280 

0.01262 
-0.00922 

1.14000 

-1.45900 

.88740 

-.40980 


St 

estimates  at  the  371    iteration,  it  can  be  seen  from  Figure 
8.12  that  the  AKI  has  essentially  converged  by  the  28 
iteration.   We  note  also  the  characteristic  convergence  of 
the  averaged  gains  to  steady  state  values  in  Figure  8.13. 
As  a  means  of  comparison  with  Perry's  results,  the  poles  and 
zeros  of  the  AKI  estimates  at  the  28  *"  ,  90    and  371    itera- 
tion are  plotted  in  Figure  8.14a  and  Figure  8.14b.   Perry's 
results  Figure  3.7  [Ref.  6:pp.  107,108]  for  his  lattice 
modeling  of  the  plant  represented  by  equation  (8.4)  are 
reproduced  for  convenience. 

As  was  noted  for  the  previous  cases,  the  instantaneous 
gains  appeared  erratic.  Figure  8.15,  whereas  the  averaged 
gains  converged  to  some  steady  state  value. 
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Figure  8.12.   AKI(4,5)  Computer  Coefficients,  Perry's  Example 
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Figure  8.13.   AKI(4,5)  Averaged  Gain  Values 
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Figure  8.15.   Instantaneous  Gain  Values,  AKI(4,5)  Model 

4 .   AKI  vs  Adaptive  Recursive  LMS  Filter 

Using  Feintuch's  proposed  algorithm  [Ref.  18]  and 
repeating  the  simulation  presented  in  the  rebuttal  to  Johnson 
and  Larimore  [Ref.  19],  a  comparison  was  made  between  the 
operation  of  the  AKI  and  the  Adaptive  Recursive  LMS  filter. 
The  Adaptive  Recursive  LMS  filter  required  "tuning"  of  the 
convergence  factors,  k   and  k  ,  to  what  appeared  to  be  the 
optimal  performance  features  of  (1)  fast  convergence  and 
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(2)  stable  estimates.   The  "tuning"  process  was  a  trial  and 
error  procedure  beginning  with  the  convergence  constants 
given  by  Feintuch, 

k    =   -4.3  X 10~^  (8.5a) 

a 

k,   =   -4.3  X  10""^  .  (8.5b) 

b 

The  Adaptive  Recursive  LMS  algorithm  was  then  implemented 
to  estimate  the  coefficients  of  the  plant  whose  transfer 
function  was , 

„,  >           0.05  -  0.40z~  ,o  r    ^ 

H(z)   =   —, -^  (8.6a) 

1  -  1.1314Z  ^    +  0.25z 

a   +  a  z"^ 
^     ^  (8.6b) 


1  -  b, z    -  bpZ 


The  best  response  obtained  by  trial  and  error  resulted  in 
the  convergence  constants, 

k    =   -4.3  X 10"^  (8.7a) 

a 

kj^   =   -14.3  X  10"^  .  (8.7b) 

2 

Femtuch  reported  that  after  8,192   iterations  the  estimates 

had  converged  on  the  coefficients  of  (8.6a)  with  a  .2096 
normalized  rms  error.   Using  the  convergence  constants  (8.7) 
convergence  was  essentially  reached  by  the  4,000    iteration 


2 

Femtuch  reports  the  resulting  estimates  for  the  coeffi- 
cients at  several  intermediate  iterations  from  8,192  to  65,536, 
however  8,192  was  the  minimum  number  of  iterations  reported. 


99 


It  was  felt  that  the  operation  of  the  adequately  tuned 
Adaptive  Recursive  LMS  filter  could  be  fairly  compared  with 
the  operation  of  the  AKI . 

Table  8.3  tabulates  the  performance  of  the  two  system 
identification  methods  and  Figures  8.16  and  8.17  graphically 
present  the  responses.   Figure  8.18  not  only  shows  at  what 
point  the  gains  of  the  AKI  reach  a  constant  value  but  also 

Table  8.3 

ACTUAL  VS  AKI  AND  ADAPTIVE  RECURSIVE  LMS  FILTER 
ESTIMATES  FOR  COEFFICIENTS  OF  EQUATION  8.6 


ADAPTIVE  RECURSIVE 
LMS  (k  =  3,990) 

0.05045 
-0.3990 

1.12800 
-0.2442 


provides  a  measure  of  confidence  that  the  unknown  system  has 
been  identified.   It  is  evident  from  Figure  8.19  that  both 
methods  seem  to  identify  the  poles  and  zeros  of  the  actual 
plant;  however,  the  poles  computed  by  the  AKI  are  closer. 
The  zeros  computed  by  the  AKI  are  approximately  .080  units 
farther  from  the  true  zero  than  is  the  Adaptive  Recursive 
LMS  filter  estimate.   All  aspects  considered,  the  AKI  takes 
considerably  fewer  iterations  to  arrive  at  its  estimate. 


ACTUAL 

AKI 

ESTIMAT 

*0 

0.05 

0.05107 

^x 

-0.40 

-0.4007 

^1 

1.1314 

1.13200 

^2 

-0.25 

-0.25090 
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Figure  8.16.   AKI(2,2)  Computed  Coefficients,  Feintuch  Example 
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Rec.  LMS  Algorithms  for  a  Plant  with  the 
Characteristic  Transfer  Function  of  Eqn .  (8.6) 
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5 .   AKI  Applied  to  PLL  Data  (Linear  Region) 

To  identify  the  coefficients  associated  with  the 
general  autoregressive-moving  average  representation  for  the 

phase  locked  loop,  the  input  noise  signal  power  (mean  square 

2 
value)  was  kept  at  (10  deg)  .   In  this  manner,  the  linear 

region  of  the  PLL  was  invoked.   The  input/output  data  was 

analyzed  by  the  AKI  resulting  in  the  estimates  shown  in 

Table  8.4.  Figures  8.20  and  8.21  show  the  respond  of  the  AKI 

Table  8.4 

ACTUAL  VS  AKI  ESTIMATES  OF  THE  ARMA 
REPRESENTATION  FOR  THE  PLL  (LINEAR  REGION) 


ACTUAL 

AKI  ESTIMATES 

(k  = 

82) 

AKI 

ESTIMATES  (k  = 

350) 

^0 

0.000000 

-0.0001046 

-0.000692 

^1 

0.020000 

0.0198500 

0.0199800 

^2 

-0.018890 

-0.018720 

-0.0188900 

^1 

1.980000 

1.980000 

1.978000 

^2 

-0.981110 

-0.981000 

-0.9790 

when  applied  to  the  identification  of  the  PLL  data.   It  was 
noted  that  though  the  AKI  correctly  identified  the  coeffi- 
cients of  the  linear  PLL,  the  AKI  gains  were  large  due  to 
the  weak  signals  (input/output  data)  incorporated  in  the 
measurement  vector,  H(k) . 

6.   AKI  Applied  to  PLL  Data  (Non-linear  Region) 

When  analyzing  any  non-linear  system  the  engineer 
must  bring  to  bear  all  his  analysis  techniques  on  the  problem, 
The  PLL  was  therefore  studied  using  classical  root  locus 
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Figure  8.20.   PLL ,  AKI(2,3)  Computed  Coefficients 
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Figure  8.21.   Averaged  Gain  Values,  AKI(2,3)  PLL  Model 
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techniques  to  preview  the  possible  outcomes  when  being  modeled 
by  the  AKI.   Stability  analysis  of  a  more  complex  PLL  system 
using  the  root  locus  technique  has  been  previously  presented 
in  the  literature  [Ref .  40] . 

The  root  locus  technique  can  be  applied  to  the  PLL 
presented  in  this  thesis  by  first  assuming  that  the  sin(*) 
block  of  Figure  5.7  is  a  variable  gain,  A.   This  is  not  a 
restrictive  assumption  since  the  PLL  during  operation  generally 
tracks  small  deviations.   The  loop  gain  can  be  written  by 
inspection  as, 

L(z)   =   A(.02)(l  -  .9445z-^)z-^  ^3^3^ 

(1  -  z"-^)^ 

A(.02) (z  -  .9445) z 
(z  -1)^ 

From  equation  (8.8)  the  root  locus  of  L(z)  is  drawn  as  shown 
in  Figure  8.22.   Even  though  the  root  locus  technique  is 
generally  used  when  the  signals  in  the  system  are  considered 
deterministic,  it  is  not  surprising  that  some  of  the  results 

obtained  are  nevertheless  valid.   When  a  moderately  strong 

2  2 

input  noise  identification  signal  [E{u  (u) }  <  (25  deg)  ]  was 

used,  the  pole- zero  locations  of  the  PLL  seem  to  follow  the 

classical  root  locus  analysis.   However,  when  the  input  noise 

2 
identification  signal  power  is  increased  beyond  (25  deg) 

the  pole-zero  locations  do  not  follow  the  expected  behavior 

predicted  using  the  root  locus  method  as  can  be  seen  in 

Figure  8.22. 
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Figure  8.22.   PLL  Root  Locus  Analysis  vs  Computer  Roots 
of  the  AKI  ARMA  Model 
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The  departure  of  the  pole-zero  behavior  from  what  was 
expected  was  analyzed  by  closer  investigation  of  the  linear 
terms  of  equation  (7.8). 


y(k)   =    .02au(k-l)  -  . 01889au (k-2)  +  ( 2 . 0  -  . 02a) y (k-1) 
lin 

+  (.01889a  - l)y(k-2)  (8.9) 


Recalling  that  a  third  order  Taylor  series  approximation  of 
the  sine  is  the  functional  expressed  by  equation  (7.7), 

sin(x)   ::   ax  +  6x^  (7.7) 

one  notes  that  the  linear  region  is  described  when  a  =  1  and 
3=0.  It  is  therefore  reasonable  to  study  the  variation  of 
a  with  respect  to  input  noise  power. 

The  average  value  a  of  a,  was  computed  by  equating 
the  estimated  AKI  coefficients  to  the  coefficients  of  like 
terms  in  equation  (8.9)  for  several  input  noise  power  levels. 
The  relation  used  was 

_     ,  2.0  -b,    1.0  +hy  a,      a~ 

^      ^      4^   :"02    "^  .01889   "^  702"  "  .01889^  (8.10) 

at  the  350th  iteration.   The  relationship  between  a  and  the 
input  noise  power  level  is  readily  apparent  from  Figure  8.23. 
This  result  is  plausible  since  if  one  considers  the  input- 
output  relationship  of  the  sine  block  of  Figure  5.7  one  obtains 
Figure  8.24.   Superimposing  the  gaussian  probability  functions 
of  the  different  input  noise  signals,  it  can  be  seen  that  for 
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low  power  levels  of  the  input  signal,  identification  of  the 

linear  parameters  of  the  overall  system  (equation  5.9)  can 

2 
be  made.   As  the  input  power  is  increased  beyond  (25  deg) 

the  input-output  characteristics  of  the  sin  x  block  are  no 
longer  approximately  linear,  showing  its  effect  in  Figure  8.24 
as  a  departure  from  its  linear  operation,  a  =  1.   Therefore, 
by  monitoring  a  one  can  determine  when  the  overall  system  is 
entering  its  non-linear  operating  regime. 

Since  the  same  functional  dependence  between  each  of 
the  AKI  estimates  and  a  did  not  exist  for  all  of  the  input 
noise  powers  considered,  knowing  a  did  not  provide  any  infor- 
mation of  the  pole  locations  but  did  provide  a  measure  of  the 
degree  of  non-linear  operation. 

B.   ORDER  OF  THE  UNKNOWN  SYSTEM  IS  NOT  KNOWN  (OVERMODELING) 
1.   AKI  vs  LMS  Adaptive  Filter 

Using  the  data  derived  by  operating  a  plant  with  the 
transfer  characteristic  of  equation  (8.1)  the  operation  of 
the  AKI  and  the  LMS  adaptive  filter  was  compared  when  the 
orders  used  in  the  identifier  and  the  LMS  filter  were  greater 
than  the  known  process  that  generated  it.   The  following 
overmodeling  cases  were  studied: 

(1)  MA  model  is  greater  than  plant  MA  process 

(2)  ARMA  model  is  fitted  to  MA  plant 

When  the  MA  model  order  is  greater  than  that  of  the 
plant  MA  process,  it  was  found  that  both  the  AKI  and  the  LMS 
filter  would  compute  coefficients  close  to  zero  for  the  higher 
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order  coefficients.   This  effect  can  be  seen  in  Figures  8.25a 
and  8.25b  when  a  fourth  order  MA  process  is  used  to  model  the 
actual  second  order  process  represented  by  equation  (8.1) . 
Since  these  are  a  "zeros  only"  plant  and  model  the  over- 
modeling  essentially  causes  zeros  to  appear  at  the  origin  of 
the  Z-plane  of  the  model  transfer  function.   Table  8.5  sum- 
marizes the  results  for  two  overmodeling  conditions  at  the  250th 
iteration  using  the  AKI .   When  a  purely  moving  average  process 
represented  by  equation  (8.1)  was  modeled  as  an  ARMA  process, 
one  must  direct  his  attention  to  the  poles  and  zeros  of  the 
model  transfer  function  which  the  AKI  computed  and  compare 
them  to  the  actual  plant  poles  and  zeros.   It  was  not  readily 
apparent  from  the  resulting  coefficients  that  the  plant  had 
been  identified.   The  following  example  will  help  clarify 
what  is  happening. 

For  the  data  produced  by  the  second  order  plant 
(equation  8.1),  an  autoregressive  moving  average  (ARMA)  model 
of  orders  2  and  5  respectively  was  fitted.   It  can  be  seen 
from  Table  8.6  that  no  firm  conclusions  can  be  drawn  about  which 
coefficients  actually  identify  the  plant.   We,  the  analysts, 
knowing  the  form  of  the  plant  which  produced  the  data,  could 
qualitatively  state  that  b,  ,  b_,  a-,  a   and  aj.  are  small 
enough  to  be  ignored.   Hence,  we  can  identify  the  plant  cor- 
rectly.  However,  inspection  of  the  poles  and  zeros  of  the 
transfer  function  which  the  AKI  computed. 


111 


3. 


UJ 
Q 


U 
11- 
li_ 
UJ 

s  • 

I 

3" 


0.       - 


-I — I 1 — I ' 1 — I 1 — I 1 1 — I 1 — 1 1 1 — I i — I 1 1 — l-H — (— I — (— I 1 — I 1 1 1 1 1 — ' 1— < — I — I i K 


oaoooogiiS'S'sooijoiDOaiiio^'Sffl'BOO^P'^'OoooooO'^aoooo 


"  0.  10. 


'''■''' 


20. 


30. 


40. 


50. 


Figure  8.25a.   AKI(0,4)  Model  of  a  Third  Order  MA  Process 
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Table  8.5 
OVERMODELING  OF  A  MA  PROCESS  USING  THE  AKI  (k  =  250) 


Actual  Model 

3rd  Order  AKI 

5th  Order  AKI 

Coeffs 

MA  Model 

MA  Model 

^0 

1.0 

.9995 

.9991 

^ 

2.0 

2.0000 

2.000 

^2 

3.0 

3.0000 

3.000 

^ 

0.0 

.893  xio"^ 

.0003191 

^4 

0.0 

.001082 

^5 

0.0 



.0003198 

zeros: 

-1.0  ±1.414j 

-1.000  ±1.414j 

-1.001  ±1.414j 

-.2977  X  10"^ 

.0227  ±  .0428J 
-.4534 

Table  8.6 

OVERMODELING  OF  A  MA  PROCESS  USING  THE  AKI  OF  ORDERS 

AR  =  2,  MJ!^  =  6 


Actual  Plant  Coefficients 

AKI(2, 

,6)   Model  Coefficients 

'^l 

0.0 

.00467 

^2 

0.0 

.10500 

^ 

1.0 

1.0000 

h 

2.0 

1.9940 

^ 

3.0 

2.8840 

^ 

0.0 

-0.2257 

^ 

0.0 

-0.3142 

^ 

0.0 

0.0001517 
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1.0    +    1.994Z"-'-   +    2.8840z"^    -    0.2257z"^ 


-    0.3142z"'^    +    0.001517z"^ 
H^T^T  ,n    ^x  (z) 


^^^2,6)  3_    _    ,00467z"^    -    .1050Z-2 


results    in: 


(8.11) 


poles:       -.8217 

.3264 
zeros:       -.9996    ±    1.414J 

.4827  X 10"^ 
-.3212 

.3262 

The    observation   to  be  made    is    that   there    are    two  pole -zero 
combinations   which    are   near   cancellation.      This    suggests    that 
the   AKI    algorithm  be    rerun  with    the   AR  and   MA  orders    reduced 
by   at   least    two. 

The    implication  of    the    analysis    of   this    section  is 
that   a  model    can   in    theory   be    found    for    a   given   set   of    input/ 
output  data.      Further,    a  parsimonious   model   can  be    identified 
by   careful   observation   of  pole-zero   combinations   which   are 
near   cancellation    and   of   stray    zeros    near    the    origin. 
2.       AKI    Applied    to   AR   and   ARMA  Data 

Essentially    the   same    characteristic   results    found   in 
Section  VIII. B.l  were    confirmed  when   the   data  produced  by 
plants   defined   by    the    transfer    functions    of   equations     (8.3) 
and    (8.4)    were    analyzed  by    the    overmodeled  AKI.      That   is, 
pole-zero  pairs    near   cancellation   and    zeros   near   the    origin 
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were  produced  by  the  AKI .   Figures  8.26  through  8.28  have 
been  chosen  as  representative  pole-zero  plots  of  the  transfer 
functions  computed  by  the  AKI  for  the  AR  plant  of  Section 
VIII. A. 2  and  for  Perry's  model  (Section  VIII. A. 3). 

C.   CONCLUSION 

This  work  indicates  that  the  Kalman  filter  algorithm 
heretofore  generally  used  as  a  state  estimator,  or  in  augmented 
form  to  estimate  parameters  (in  which  case  the  parameters 
are  treated  as  states) ,  can  also  be  formulated  in  an  adaptive 
manner  to  iteratively  estimate  the  coefficients  of  an  ARMA 
equation  explicitly.   This  approach,  termed  the  Adaptive 
Kalman  Identifier  (AKI) ,  summarily  identifies  the  unknown 
system  whose  input/output  data  is  being  processed.   The  LMS 
adaptive  algorithm  of  Widrow,  and  its  modification  by  Griffiths 
(in  which  the  convergence  factor  is  selected  to  be  inversely 
proportional  to  the  input  signal  power)  are  shown  to  be  sub- 
optimal  cases  of  the  AKI.   An  additional  insight  provided  by 
the  AKI  is  that  it  indicates  clearly  how  measurem.ent  noise 
might  be  taken  into  account  in  the  LMS  adaptive  formulation. 

The  operation  of  the  AKI  was  checked  by  way  of  simulation 
and  compared  with  two  existing  identification  techniques: 
(1)  the  LMS  Adaptive  nonrecursive  (MA)  algorithm  and  (2)  the 
Adaptive  Recursive  ARMA  LMS  algorithm.   It  was  found  that  not 
only  are  the  two  LMS  filtering  techniques  special  suboptimal 
cases  of  the  AKI;  but,  further,  the  AKI  exhibits  superior 
convergence  and  modeling  properties  for  the  cases  where  (1)  the 
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Figure    8.26.       AKI(3,1)    Overmodeled  AKI (4 , 2) 
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Figure    8.27.      AKI (3,1)    Overmodeled   AKI (5,5) 
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Figure    8.28.      AKI(4,5)    Overmodeled   AKI (5 , 6) 
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order  of  the  unknown  plant  is  known  and  (2)  the  order  of  the 
plant  is  overmodeled.   Additionally,  the  simulations  indicate 
that  accuracies  similar  to  those  obtained  using  lattice 
modeling  methods,  can  be  achieved  using  the  AKI  at  a  decidedly 
smaller  number  of  iterations. 

By  making  minor  modifications  to  the  measurement  vector, 
H(k) ,  that  is  by  using  hybrid  signals,  the  AKI  was  used  to 
identify  the  linear  and  non-linear  ARMA  representations  of  a 
phase  locked  loop  with  success.   Interestingly,  the  AKI  tech- 
nique appears  to  enable  one  to  discern  when  a  potential  non- 
linear system  enters  its  non-linear  mode  of  operation,  by 
closely  monitoring  the  coefficients  of  the  linear  portion  of 
the  generalized  non-linear  ARMA  model. 

D.   TOPICS  FOR  FURTHER  CONSIDERATION 

Several  areas  for  further  study  directly  and  indirectly 
related  to  the  AKI  were  uncovered.   Foremost,  a  rigorous 
convergence  proof  is  desirable.   Although  the  connection  was 
made  between  the  AKI  algorithm  and  the  initial  equation  from 
which  the  lattice  modeling  algorithm  is  developed,  similar 
comparisons  (such  as  Chapters  III  and  IV)  could  provide  more 
insight  into  the  operation  of  both.   The  multichannel  AKI  is 
the  logical  development  of  the  single  channel  AKI  presented. 
And,  lastly,  refinement  of  the  A-KI  software  (using  the  NL  =  5 
option)  to  include  ARMA  modeling  of  time  series  using  a  Box- 
Jenkins  approach  [Ref.  26]  is  feasible.   A  limited  number  of 
simulations  using  Monterey  rain  data  and  series  C  [Ref.  26] 
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data  seem  to  indicate  that  an  application  exists  for  the  AKI 
in  this  area. 
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APPENDIX  A 
THE  DISCRETE  WEINER  PROBLEM 

The  discrete  Weiner  problem  considers  optimally  filtering 
a  desired  signal  from  unwanted  stationary  noise.   The  cri- 
terion of  optimality  used  is  minimization  of  the  mean  squared 
error  between  the  output  and  the  desired  signal.   Generally 
one  desires  that  the  optimal  filter,  which  is  the  device 
being  sought,  be  time  invariant  and  that  it  be  able  to  accept 
a  signal,  s(k),  and  noise,  n(k),  where  each  are  samples  from 
stationary  random  processes.   In  other  words,  we  want  a 
device  which  can  accept 

x(k)   =   s(k)  +  n(k)  (A.l) 

as  an  input  and  produce  at  its  output  s(k+A)  or  some  linear 
function  thereof  where  A  is  some  known  delay. 

Assuming  that  the  desired  output,  d(k),  is  the  response 
to  a  specified  sampled  data  linear  system  whose  transfer 
function,  H,  (z),  is  given,  then  the  error  between  the  desired 
output  and  the  output  of  the  filter,  y(k) ,  we  seek  is, 

e(k)   =  d(k)  -  y(k)  .  (A. 2) 

Figure  A.l  depicts  the  discrete  Weiner  problem  formulation. 
The  derivation  from  here  follows  the  one  presented  by  Maybeck 
[Ref.  27]  for  the  continuous  case.   From  the  block  diagram 
we  have , 
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s(k) 


n(k) 


1    h^(k-i)[-] 

]_=— CO 


d(k) 


^    £(k) 


Figure  A.l.       The    Discrete   Weiner    Problem 


Y(k)       =  I    h-(k-i)    x(i) 


(A. 3) 


or   equivalently , 


y(k)       =  I    x(k-i)    h    (i) 


(A. 4) 


Substituting  equation  A. 4  into  A. 2  and  squaring  we  have. 


e^(k)   =   d^(k)  -  2d(k)    I    h.(i)x(k-i) 

i=-oo   ^ 


(A. 5) 


+  [_  I     x(i)h^(k-i)]  [  I     x(i)h  (k-i)] 

i  =  — 00  j_=  — 00 
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Taking   the   expected  value   of   A. 5   we   obtain. 


E{e^(k)}      =      ^J^^iO)    -    2        I    h    (i)-^       (i) 


+    [  h^(i)     I  h^{k)lp^^{k-i) 

j_:=— 00  1^:=— oo 


(A. 6) 


where  notation  \b       (m)  denotes  the  expected  value  of  the 

uv  ^ 

product  of  u(k)  and  v(k+m) .   That  is. 


i)       (m)   =   E{u(k)v(k+]T>)  }   =   E{u  (k)  v  (k-m)  }         (A.  7) 


otherwise  known  as  the  autocorrelation  of  u(k)  and  v(k) . 
Using  variational  techniques  and  letting, 

h.(i)   =  h_^(i)  +  £Ah(i)  (A. 8) 

f         opt 

substituting  equation  A. 8  into  equation  A. 6,  taking  the 
partial  derivative  with  respect  to  epsilon,  e,    and  setting 
the  partial  derivative  equal  to  zero  we  have: 

OO  00 

I    Ah(n)[  I    h^^.(i)i|;  ^(n-i)  -  4^  .(n)]   =   0        (A. 9) 
^        .  ^   opt    XX         xd 


n=-oo 


.  ^   op 

1=  — OO     ^ 


The  term  within  the  brackets  is  the  discrete  form  of  the 
Weiner-Hopf  equation  and  must  be  equal  to  zero  if  equation  A. 9 
is  to  be  valid  since  Ah(n)  >  0  by  definition.   Therefore, 

00 

I    l^^^^(i)'4^   (n-i)   =   i|^  .(n)  (A. 10) 

.  ^   opt    XX  xd 

J_^  — 00    '- 
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Equation  A. 10  is  the  most  often  encountered  in  truncated  form 
in  linear  prediction  theory.  Fade  approximation,  adaptive 
filtering  and  lattice  filtering.   The  truncated  version  of 
A. 10  which  is  generally  used  in  solving  the  discrete  Weiner 
problem  is. 


M 


i  =  0 


h   ^(i)i|;   (n-i) 

opt      XX 


=   '^xd^^) 


(A. 11) 


or   in  matrix   form. 


ip       (0)       ij;       (-1) 

XX  XX 


^       (-n) 

XX 


i>       (1)       HJ       (0)  ...       ^i)       (1-n) 

^XX  XX  ^xx 


\l)       (n)       ^      (n-1)     ...      ip       (0) 
^xx  ^xx  ^xx 


h      ^(0) 
opt 

opt 


J    L  opt 


*xd("' 


(A. 12) 
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